Modelling the Reverse ElectroDialysis process with seawater and concentrated brines

This article was downloaded by: [Universita di Palermo], [Andrea Cipollina]On: 27 November 2012, At: 10:00Publisher: Taylor & FrancisInforma Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House,37-41 Mortimer Street, London W1T 3JH, UK

Desalination and Water TreatmentPublication details, including instructions for authors and subscription information:http://www.tandfonline.com/loi/tdwt20

Modelling the Reverse ElectroDialysis process withseawater and concentrated brinesMichele Tedesco a , Andrea Cipollina a , Alessandro Tamburini a , Willem van Baak b &Giorgio Micale aa Dipartimento di Ingegneria Chimica, Gestionale, Informatica, Meccanica, università diPalermo, viale delle Scienze (Ed.6), 90128, Palermo, Italyb FUJIFILM Manufacturing Europe BV, Oudenstaart 1, P.O. Box 90156, 5000LJ, Tilburg, TheNetherlandsVersion of record first published: 10 Jul 2012.

To cite this article: Michele Tedesco , Andrea Cipollina , Alessandro Tamburini , Willem van Baak & Giorgio Micale (2012):Modelling the Reverse ElectroDialysis process with seawater and concentrated brines, Desalination and Water Treatment,49:1-3, 404-424

To link to this article: http://dx.doi.org/10.1080/19443994.2012.699355

PLEASE SCROLL DOWN FOR ARTICLE

Full terms and conditions of use: http://www.tandfonline.com/page/terms-and-conditions

This article may be used for research, teaching, and private study purposes. Any substantial or systematicreproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form toanyone is expressly forbidden.

The publisher does not give any warranty express or implied or make any representation that the contentswill be complete or accurate or up to date. The accuracy of any instructions, formulae, and drug doses shouldbe independently verified with primary sources. The publisher shall not be liable for any loss, actions, claims,proceedings, demand, or costs or damages whatsoever or howsoever caused arising directly or indirectly inconnection with or arising out of the use of this material.

http://www.tandfonline.com/loi/tdwt20

http://dx.doi.org/10.1080/19443994.2012.699355

http://www.tandfonline.com/page/terms-and-conditions

Modelling the Reverse ElectroDialysis process with seawater andconcentrated brines

Michele Tedescoa, Andrea Cipollinaa,*, Alessandro Tamburinia, Willem van Baakb,Giorgio Micalea

aDipartimento di Ingegneria Chimica, Gestionale, Informatica, Meccanica, universita di Palermo, viale delle Scienze(Ed.6), 90128 Palermo, ItalyEmail: [email protected] Manufacturing Europe BV, Oudenstaart 1, P.O. Box 90156, 5000LJ Tilburg, The Netherlands

Received 12 March 2012; Accepted 29 May 2012

ABSTRACT

Technologies for the exploitation of renewable energies have been dramatically increasing innumber, complexity and type of source adopted. Among the others, the use of saline gradi-ent power is one of the latest emerging possibilities, related to the use of the osmotic/chemi-cal potential energy of concentrated saline solutions. Nowadays, the fate of this renewableenergy source is intrinsically linked to the development of the pressure retarded osmosisand reverse electrodialysis technologies. In the latter, the different concentrations of two sal-ine solutions is used as a driving force for the direct production of electricity within a stackvery similar to the conventional electrodialysis ones. In the present work, carried out in theEU-FP7 funded REAPower project, a multi-scale mathematical model for the Salinity Gradi-ent Power Reverse Electrodialysis (SGP-RE) process with seawater and concentrated brineshas been developed. The model is based on mass balance and constitutive equations col-lected from relevant scientific literature for the simulation of the process under extreme con-ditions of solutions concentration. A multi-scale structure allows the simulation of the singlecell pair and the entire SGP-RE stack. The first can be seen as the elementary repeating unitconstituted by cationic and anionic membrane and the relevant two channels where diluteand concentrate streams flow. The reverse electro-dialysis stack is constituted by a numberof cell pairs, the electrode compartments and the feed streams distribution system. Themodel has been implemented using gPROMS�, a powerful dynamic modelling process sim-ulator. Experimental information, collected from the FUJIFILM laboratories in Tilburg (theNetherlands), has been used to perform the tuning of model formulation and eventually tovalidate model predictions under different operating conditions. Finally, the model has beenused to simulate different possible scenarios and perform a preliminary analysis of the influ-ence of some process operating conditions on the final stack performance.

Keywords: Saline gradient power; Reverse electrodialysis; Modelling; Multi-scale; gPROMS

*Corresponding author.

Presented at the International Conference on Desalination for the Environment, Clean Water and Energy, European DesalinationSociety, 23–26 April 2012, Barcelona, Spain

Desalination and Water Treatmentwww.deswater.com

1944-3994/1944-3986 � 2012 Desalination Publications. All rights reserveddoi: 10.1080/19443994.2012.699355

49 (2012) 404–424

November

Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

1. Introduction

The increasing world energy demand during lastdecades, in conjunction with sustainability issuesrelated to the large use of fossil fuels, is rapidly lead-ing to a significant interest towards research on newpossible renewable energy sources. Among the severalinvestigated researches, a promising option is the useof salinity gradient power (SGP), also called “osmoticenergy”, i.e. the chemical energy potential associatedwith the “controlled” mixing of two salt solutions atdifferent concentrations.

The concept of SGP is already known in the litera-ture and was described for the first time in 1954 byPattle [1]. Two different techniques have been pro-posed in the literature for recovering the osmoticenergy of a system and convert it into a more exploit-able form: pressure retarded osmosis (SGP-PRO) andreverse electrodialysis (SGP-RE). In the former, SGP isconverted into mechanical energy (and then to electri-cal energy by means of turbines) using osmotic mem-branes, i.e. membranes allowing the passage of waterand obstructing the passage of salts; in the latter, SGPis converted directly into electric energy adoptingselective ion exchange membranes (IEM) within astack similar to conventional electrodialysis ones.

In general, the overall process performance isstrongly related to the concentration gradient avail-able, e.g. river water, seawater or differently concen-trated saline solutions. In particular, recently, therehas been a growing interest for the use of salinity gra-dients generated by the mixing of seawater and veryconcentrated brines such as those produced in salt-works, salt mines or some other industrial activities.In this regard, it is worth mentioning the EU-FP7funded REAPower project [2], which aims at thedevelopment of a SGP-RE unit for power productionfrom seawater and concentrated brines.

1.1 Principle of SGP-RE

A simplified scheme of a typical SGP-RE stack isshown in Fig. 1. The key components of the system arethe ion exchange membranes (IEMs), which are assem-bled in a stack with alternating cationic-type membrane(CEM) and anionic-type membrane (AEM). The dis-tance between two subsequent membranes (i.e. thecompartment thickness) is generally guaranteed by theuse of polymeric spacers.

In particular, the repeating unit of the stack (gener-ally called cell pair) consists of four elements: (i) aCEM; (ii) a compartment for the flowing concentratedsolution (concentrate compartment); (iii) an AEM and (iv) acompartment for the flowing diluted solution (diluate

compartment). The two solutions are forced to flowthrough the stack in alternate channels by means ofsuitably shaped inlet and outlet distribution systems.The external compartments of the stack (electrode com-partments) contain the electrodes and the flowing electro-lyte solution (electrode rinse solution), which also containsa suitable redox couple (e.g. Fe2+/Fe3+ chloride).

In the scheme shown in Fig. 1, the two solutionsare fed to the stack in co-current mode: in general, inmembrane separation process (such as conventionalelectrodialysis process) it is preferable to use this con-figuration just to avoid significant local pressure dif-ference between the two sides of IEMs and preventthe risk of internal leakages and mechanical stress tomembranes [3].

With reference to the scheme reported in Fig. 1,the principle of SGP-RE process can be synthesised inthe following: when the two salt solutions are fed tothe stack, the concentration gradient between concen-trate and diluate acts as driving force for the ions todiffuse across membranes. The passage of ionsthrough a membrane is regulated by its permselectivi-ty, i.e. the selectivity towards the passage of positivelyor negatively charged ions rather than others (e.g. cat-ions pass through CEM, while anions are rejected). Inideal conditions, cations (mainly Na+ ions in seawatersolutions) flow only across CEMs, while anions(mainly Cl� ions) pass only through AEMs, thus mov-ing in opposite directions. These ionic fluxes acrossmembranes constitute the ionic current through thestack, which is eventually converted into electric cur-rent at the electrodes. In fact, the role of the electroderinse solution is to restore electroneutrality in externalchannels by means of redox reactions at the elec-trodes: in this way the electric continuity of the sys-tem is ensured and the generated electric current canbe used by an external load.

1.2 State–of-the-art of modelling SGP-RE process

The earliest works regarding theoretical calcula-tions on the SGP-RE process have been proposed byWeinstein and Leitz [4] and Forgacs and O’Brien [5].However, these works, which also represent the firstexperimental demonstration of SGP-RE process techni-cal feasibility, only presented the use of standard ther-modynamics and “electrical model” equations for theestimation of theoretical obtainable power output.Moreover, results are certainly affected by the use ofimportant simplifying assumptions relevant tosolutions/membranes properties (e.g. the use ofconcentrations rather than activities, or neglecting theeffects of concentration on physical properties).

M. Tedesco et al. / Desalination and Water Treatment 49 (2012) 404–424 405

Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

The first comprehensive work on SGP-RE proc-ess modelling was proposed by Lacey [6], whosuggested a novel approach for estimating the cellpair voltage as the sum of electrical potential dropsbetween the two bulk solutions through the polari-sation layers, diffusive boundary layers and IEMs.This work, though lacking of any correlation forphysical properties calculation, clearly demonstratedthe importance of a careful equipment design forimproving process efficiency. With this regard,Lacey underlined the importance of using a largeconcentration ratio between concentrate/dilute solu-tions and suitably selecting diluate concentration inorder to prevent high electric resistance in dilutecompartment.

Lacey’s model was later modified by Brauns [7]and implemented in a solver software in order toinvestigate the effect of specific parameters on processefficiency. In particular, Brauns analysed thepossibility of using two different thicknesses fordiluate/concentrate compartments, investigating

especially the effect of reducing electrical resistance byusing a thinner diluate compartment. Furthermore,the effects of membrane thickness, inlet solutions con-centration and temperature on cell pair voltage wereinvestigated. Simulation results show how power out-put is largely affected by membrane thickness, indicat-ing that the development of novel thinner membranescan give a significant enhancement of processefficiency.

Finally, a model for SGP-RE cell pair was pro-posed by Veerman et al. [8], based on thermodynamicequations for cell pair voltage calculation and trans-port equations coupled with electrical model equa-tions for mass/charge flux across IEMs. The modelalso included three adjustable parameters, physicallyrelated to spacer shadow effects, osmotic diffusionand co-ionic flux through the membrane. A completeset of mass balance equations were also included inthe model to take into account also for the variation ofsalt concentration along the channel due to co- andcounter-ions fluxes through the membranes.

AEM CEM AEMCEM AEMCEM

Electrode Rinse Solution

Electrode Rinse Solution

CEM

BRINE INSEAWATER IN

BRINE OUT

Cell Pair

e-e-

External load

CATHODE(+)

ANODE(-)

+Na+

-Cl -

+Na+

-Cl-

+Na+

-Cl-

Red

Ox

Red

Ox

SEAWATER OUT

Fig. 1. Simplified scheme of a SGP-RE unit.

406 M. Tedesco et al. / Desalination and Water Treatment 49 (2012) 404–424

Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

Two different “efficiency parameters” were takeninto account for investigating the effect of varyingoperating conditions: (1) the “net power density”,defined as the obtainable net power (power outputminus pumping power required for pumping inletsolutions) normalised with respect to total membranearea and (2) “riverwater yield”, defined as the ratiobetween net power and required diluate flow rate.This latter is defined assuming that the availability ofriver water can be a limiting factor for the process.After a suitable tuning of model adjustable parametersby comparison of purposely collected experimentalinformation, an optimisation study was carried outwith respect to these output variables, analysing theeffects of compartments thickness and residence timeof solutions.

On the basis of the literature review it seems clearhow, although some efforts have been done alreadyfor implementing models for the SGP-RE process,these works are only related to the use of river waterand seawater as feed solutions, while none of themhas considered the use of concentrated brines, whichcould dramatically affect the model formulation.Moreover, all presented approaches are based on anumber of simplifying assumptions, which could beovertaken for example by the use of a “multi-scale”approach, where different modelling scales are usedto simulate different aspects of the process such asfluid flow behaviour, heterogeneity of flows distribu-tion, transport phenomena inside cells and electricalphenomena within the entire stack.

1.3 Focus of the work

Focus of this paper is to present a multiscalemodel approach for the simulation of a SGP-RE stackusing “seawater” and “concentrated brine” as feedsolutions. In particular, an aqueous solution with aconcentration of 30 g/l of NaCl (0.52M) was consid-ered as diluate and a 315 g/l NaCl aqueous solution(5.4M) as concentrate stream. This last value isdirectly related to the real salinity of concentratedbrine let out from Trapani’s saltworks, where a proto-type will be installed within the activities of theREAPower (Reverse Electrodialysis Alternative Power)EU-funded project [2].

The use of a multi-scale modelling approach isparticularly suitable due to the noteworthy complexityof the process and the major goals of the modellingitself: i.e. the development of an advanced modellingtool for the optimisation of the process and the designof a SGP-RE stack to be constructed and tested withinthe REAPower project.

Indeed, the multiscale approach allows the mathe-matical description of the system at different scales, fromthe single cell pair to the entire stack and, eventually, tothe complete plant including also auxiliary units.

Following the analysis of the state-of-the-art onSGP-RE process modelling, the cell pair model pro-posed in 2011 by Veerman et al. [8] for river water/seawater SGP-RE has been considered and re-adaptedto the more complex case of seawater and brine. Then,starting from the cell pair model, a hierarchical struc-ture has been implemented using gPROMS, a power-ful dynamic process simulator, aiming at themodelling of the entire stack behaviour, which hasbeen eventually validated by comparison with experi-mental information.

2. Model development

2.1. Development of the lower hierarchy model (cell pair)

2.1.1. System definition

A sketch of a single cell pair in co-current modeoperation is shown in Fig. 2. System geometry can becharacterised by four main parameters: membranelength (L) and width (b), spacer thicknesses of diluatecompartment (ds) and concentrate compartment (db).

Simplifying assumptions have been done formodel implementation [8]:

(a) both solutions are assumed to be constitutedonly by sodium chloride (i.e. neglecting all ofother ions present in real seawater/brine solu-tions);

(b) salt concentration profiles (in the direction per-pendicular to the IEMs) are assumed to be flatin both compartments: i.e. concentration polari-sation phenomena are neglected;

(c) permselectivity of both IEMs is assumed con-stant;

(d) electro-osmotic flux is considered negligible;(e) a bi-dimensional approach has been used, thus

neglecting any possible variation in the direc-tion of channel width.

In order to take into account the non-ideality ofmembranes two “disturbing effects” are considered:(1) co-ions transport through the membranes and (2)solvent osmotic flux through the membranes.

The former is related to the diffusion of ions hav-ing the same charge of membrane functional groups(e.g. cations for AEMs) driven only by the concentra-tion difference between the membrane sides. In thisway, co-ions move against the main current direction,


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

thus having a negative effect on process efficiency. Onthe other side, the effect of osmotic flux (from diluateto concentrate compartment) leads to a furtherreduction of concentration gradient between concen-trate–diluate and, thus, of driving force for powerproduction.

A “distributed parameters” model has been imple-mented, considering the variation of main variables(concentrations, activities, cell potentials, etc.) alongthe channel length (L), where a x-coordinate has beendefined. Finally, all specific variables refer to the cellpair membrane area.

2.1.2. Model equations

Constitutive equations for thermodynamic properties.Veerman et al. [8] estimated activity coefficients bythe Extended Debye–Huckel (EDH) equation. Thiswell-known theory is rigorously valid only for electro-lyte solutions at infinite dilution. In particular, for 1:1-valent electrolytes (such as NaCl) in aqueous solutionthe EDH theory gives a very good accuracy in predict-ing experimental trends up to a normal concentrationof 0.5 eq/L [9], which is comparable with the normalconcentration of seawater.

In order to estimate mean activity coefficients forelectrolyte in more concentrated solutions, where con-centrations range from 0.5 eq/L up to values close to

saturation, a virial equation proposed by Pitzer [10]has been used:

ln c� ¼ fc þmBc þm2Cc ð1Þ

where c� and m are the mean activity coefficient andmolality of the electrolyte, respectively. The term (f c),taking into account long-range forces between ions, isdirectly derived from EDH theory; while the secondand third virial coefficients (Bc, Cc) are related toshort-range interactions. For the estimation of thesetwo binary interactions parameters (a, b) are adopted,which depend on the nature of the electrolyte and canbe easily found in the open literature for the mostcommon substances. For the sake of brevity, relevantexpressions are not reported, but can be easily foundthe literature [10]. To confirm what is stated above,Fig. 3 shows a comparison between predicted andexperimental values of mean activity coefficient forNaCl as a function of electrolyte concentration. It isworth noting how, for molar concentration above0.5M, EDH theory, predicting a monotone trend forc�, is not able to give accurate values of activity coeffi-cient, while a good agreement with experimentaltrend is observed for the Pitzer equation.

Another physical property to be estimated formodelling purposes is the equivalent conductivity. Itcan significantly vary with concentration, especiallywhen it is far from the infinite dilution condition. Forthis reason, while in Ref. [8] equivalent conductivitywas assumed constant, in the present work concentra-tion influence on conductivity has been considered.

Among the various correlations proposed in litera-ture, the most accurate for NaCl solutions is the onesuggested by Islam et al. [12,13], which is based onthe Debye–Huckel–Onsager theory:

� ¼ �0 � B02ðcÞ

ffiffic

p1þ B0ðcÞa ffiffi

cp

� �1� B0

1ðcÞffiffic

p1þ B0ðcÞa ffiffi

cp F0ðcÞ

� �ð2Þ

where K0 is the equivalent conductivity at infinitedilution and c is electrolyte molar concentration; theother terms of Eq. (2) (whose definitions are notreported herein for the sake of brevity) are functionsof concentration as well as some physical propertiesof solution itself, such as viscosity (g) and dielectricconstant (e). Thus, in order to evaluate K reliable datafor g and e at different salt concentrations are needed:e.g. for a 5M NaCl solution e is roughly 50% lowerthan for pure water [9]. In this paper, the effect ofelectrolyte concentration on g and e has been consid-ered assuming a linear trend with salt concentrationfor both [11].

Fig. 2. Simplified scheme of a repeating unit within a SGP-RE stack (cell pair).


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

Fig. 4 shows the estimated equivalent conductivityfor NaCl solution at different molar concentration;note that, assuming a constant value for both g and e,the accordance with experimental data is acceptableonly for concentration lower than 1M, while a verygood matching can be achieved with the assumedlinear dependence of these parameters with the con-centration.

Electric variables. The open circuit voltage (OCV) ofthe system, i.e. the potential difference created at bothIEMs, can be calculated by an equation formallysimilar to the Nernst equation [8]:

EcellðxÞ ¼ aCEMRT

Fln

cNab ðxÞCbðxÞcNas ðxÞCsðxÞ þ aAEM

RT

F

� lncClb ðxÞCbðxÞcCls ðxÞCsðxÞ ð3Þ

where aAEM, aCEM are the permselectivity of bothIEMs, and the subscripts b and s are referred to brineand seawater, respectively.

The total cell pair resistance is constituted by thesum of 4 resistances in series:

RcellðxÞ ¼ RbðxÞ þ RsðxÞ þ RCEM þ RAEM ð4Þ

where Rb and Rs are electrical resistances of solutions,while RCEM and RAEM are the IEMs resistances; in Eq.

(4) all terms are expressed as areal resistances (Xm2),as normally done for this kind of systems. As it con-cerns solutions resistances, they are calculated by:

Rb ¼ fdb

�bCbðxÞ ð5Þ

Rs ¼ fds

�sCsðxÞ ð6Þ

where f is called obstruction factor and represents a cor-rection term which takes into account the increase ofelectrical resistance caused by the presence of thespacer (e.g. due to tortuosity of ion path or shadoweffect on membranes).

Transport equations across IEMs. The concentrationgradient across each membrane generates the trans-port of ions from concentrate to dilute compartment.In general, because of the non-ideality of membranes,both counter-ions and co-ions will pass through eachIEM; as a consequence, the total salt flux across IEMscan be divided in two terms (Eq. (7)):

JtotðxÞ ¼ JcoulðxÞ þ JcitðxÞ ð7Þ

in Eq. (7) the first term (counter-ions or “coulombic”flux, Jcoul) represents the ions flux responsible of the

0

25

50

75

100

125

150

0 1 2 3 4 5 6

Equi

vale

nt c

ondu

ctiv

ity, Λ

[S c

m2 / m

ol]

NaCl molar concentration [mol/l]

experimental data

Islam equation

modified Islam equation

Λ0

Fig. 4. Equivalent conductivity of NaCl aqueous solution.Comparison between values predicted by Islam equation,modified Islam equation and experimental data.Experimental data from Ref. [11].

0.5

0.6

0.7

0.8

0.9

1.0

0 1 2 3 4 5 6

mea

n ac

tivity

coe

ffici

ent (

γ ±)

NaCl molal concentration [mol/Kg]

experimental dataPitzer equationEDH equation

Fig. 3. Influence of NaCl molal concentration on meanactivity coefficient.Comparison between values predicted by EDH, Pitzerequation and experimental data. Experimental data fromRef. [11].


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

electric current in external circuit: this contributioncan be related to the current density (j) simply by theFaraday’s law. The second term, related to co-ionstransport (Jcit), only represents a loss of driving forcefor SGP-RE process and herein is described through aphenomenological expression. Thus Eq. (7) can befurther developed as:

JtotðxÞ ¼ JcoulðxÞ þ JcitðxÞ

¼ jðxÞF

þ 2DNaCl

dm½CbðxÞ � CsðxÞ� ð8Þ

where j is the current density, dm is membrane thick-ness (assumed equal for both IEMs) and DNaCl is theco-ions diffusion coefficient; the factor 2 allows to takeinto account both IEMs of a cell pair.

In Eq. (8) the DNaCl term, formally analogous to asalt diffusion coefficient, could be considered strictlyas the mean value of co-ions diffusivities in eachmembrane (i.e. Na+ diffusivity in AEMs and Cl� dif-fusivity in CEMs, respectively). While, in Ref. [8]DNaCl was considered and adjustable parameter, inthe present model formulation, it has been estimatedfrom permselectivity definition itself. In fact, the defi-nition of CEM permselectivity commonly accepted inliterature [3] is:

aCEM ¼ TCEMc � Tc

Ta

ð9Þ

where TCEMc is the transport number of cation into the

membrane, while Tc and Ta are cation and aniontransport number in solution, respectively. Assumingonly Na+ and Cl� ions, TCEM

c should be equal to theratio between the counter-ion molar flux (Jcoul) andthe total salt flux (Jtot):

TCEMc ¼ zcJcP

i ziJi¼ Jcoul

Jcoul þ Jcitð10Þ

substituting Eqs. (8) and (10) in Eq. (9) and rearrang-ing, DNaCl can be finally expressed as:

DNaClðxÞ ¼ 1

2

1� TCEMc

TCEMc

� �jðxÞF

dmðCbðxÞ � CsðxÞÞ ð11Þ

note that, in case of ideal membrane, TCEMc ¼ 1 and

therefore co-ions diffusion coefficient is zero.

Osmotic transport. Osmotic transport of solventcan be considered as proportional to the salt concen-tration gradient between concentrate–dilute compart-ment through a phenomenological expression (Eq.(12)):

J0wðxÞ ¼ �2Dw

dm½CbðxÞ � CSðxÞ� � MH2O

qH2O

� �ð12Þ

where Dw is the water diffusion coefficient throughIEMs, MH2O and qH2O are solvent molar weight anddensity, respectively, and the multiplying factorðMH2O=qH2OÞ is added in order to express J0w as a vol-umetric flux (m3/m2 s).

Mass balance. Mass balance equations for NaCl ineach compartment can be easily expressed as:

dCsðxÞdx

¼ b

Qs

JtotðxÞ ð13Þ

dCbðxÞdx

¼ � b

Qb

JtotðxÞ ð14Þ

where Qs, Qb are the solution flow rates in singlechannel. In order to include also the effect of solventosmotic transport mass balance can be modified, asalready reported in Ref. [8]:

dCsðxÞdx

¼ b

Qs

JtotðxÞ � CsðxÞ b

Qs

J0wðxÞ ð15Þ

dCbðxÞdx

¼ � b

Qb

JtotðxÞ þ CbðxÞ b

QbJ0wðxÞ ð16Þ

indicating a further reduction of concentration in theconcentrate compartment and an increase in thediluate one, as physically expected.

2.2. Definition of the adjustable parameter

In the model proposed by Veerman et al. [8],the obstruction factor (f) and diffusion coefficientsof water and salt (Dw, DNacl) were used as adjust-able parameters. In particular, f appears within cal-culation of areal resistances of solutions (Eqs. (5)and (6)). The influence of f on overall process effi-ciency, together with its possible use as adjustableparameter for the model, could be an importantpoint to discuss. In fact when river water is usedas dilute, the areal resistance of dilute solution is


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

the controlling resistance for the system, as alreadyunderlined by several authors [7,14]; thus it is clearthat the obstruction factor has a great influence onfitting in that case.

On the other hand, concerning the conditionsinvestigated in this work (SGP-RE process with sea-water–brine) the higher salt concentration within thechannels makes the areal resistance of both solutions(Rs, Rb) negligible with respect to IEMs resistances(RCEM, RAEM). Therefore, the effect of the obstructionfactor on the overall process efficiency becomes poorand its value has been set as done by Veerman et al.(i.e. f= 2.5), after assessing its negligible influence onthe process performance. For similar reasons, also Dw

has been fixed to the literature value (i.e.Dw=1� 10�9m2/s) [8]. While DNaCl has been elimi-nated by as indicated above (Eq. (11)).

A new adjustable parameter (more suitable to thefitting) has been added to the model, in order to cor-rect nominal values of permselectivities of IEMs totake into account of the effects of high solution con-centrations, thus estimating a “corrected” OCV fromEq. (3):

EcellðxÞ ¼ b � ðaCEM þ aAEMÞRTF

lncbðxÞCbðxÞcsðxÞCsðxÞ ð17Þ

where b is the permselectivity correction factor. Notethat in this simplified case (one single electrolyte insolution), the two terms on second member in Eq. (3)can be lumped together. Notwithstanding its defini-tion, b can be seen also as a fitting parameter for tak-ing into account any non-ideal effect related to themembranes, e.g. the shadow effect generated by thespacer filaments in contact with the membranesurface.

2.3. Development of the higher hierarchy model (stack)

Once the model for cell pair has been developed,the next step is the construction of a model for theentire stack, i.e. a system constituted by a genericnumber N of cell pairs plus the electrode compart-ments. While transport phenomena through IEMs areconsidered within each cell pair model, the main vari-ables within the “stack model” are the electric vari-ables that can be estimated only after a completemathematical description of the entire system.

The total stack resistance is the sum of cells resis-tances plus the resistance in electrode compartments(blank resistance, Rblank):

RstackðxÞ ¼XNi¼1

Rcell;iðxÞ !

þ Rblank ð18Þ

of course, the contribution of electrode compartmentson total stack resistance is relevant only for stack withsmall number of cells (e.g. in lab-scale units).

The potential difference available at the stack willbe the sum of the OCV in all cells minus the potentialdrop due to the internal stack resistance:

EstackðxÞ ¼XNi¼1

Ecell;iðxÞ !

� RstackðxÞ � jðxÞ ð19Þ

Then, the current density can be calculated by Ohm’slaw:

jðxÞ ¼ EstackðxÞRu

ð20Þ

where Ru is the external load of the system (assumedas constant), expressed in Xm2.

Aside from the OCV, another useful quantity justto express the efficiency of a SGP-RE system is repre-sented by the specific power, or power density, whichis defined as the obtainable electric power per mem-brane area of cell pair. This variable represents theavailable output power normalised with respect to therequired membrane surface (which is one of the prin-cipal cost issue for the process itself) and can be calcu-lated as:

PdðxÞ ¼ 1

Nj2ðxÞ � Ru ð21Þ

note that in Eq. (21) power density is defined withrespect to the cell pair area.

The output electrical power obtainable from thestack (Pout), is then calculated by Joule’s law as prod-uct of average stack voltage and average electric cur-rent (Pout ¼ �Estack � �I) or by multiplying Pd by the totalstack cell pair area.

The total required power for pumping the feedsolutions through the stack can be estimated as:

Ppump ¼ DPb �Qtotb þ DPs �Qtot

s

gpð22Þ

where DP is the total pressure drop (Pa), Qtot is thesolution flowrate (m3/s), gp is the pump efficiency,subscripts b and s refer to brine and seawater,respectively.


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

Hence, the net power, which is the obtainable elec-tric power subtracting the required pumping power,can be estimated as:

Pnet ¼ Pout � Ppump

¼ ð�Estack � �IÞ � DPb �Qtotb þ DPs �Qtot

s

gpð23Þ

2.4. Model implementation and simulating approach

The proposed DAEs set has been implemented ingPROMS�, a powerful software for simulation (bothin steady state and dynamic conditions) and optimisa-tion of complex processes [15]. Furthermore, unlikecommon flow-sheeting software, which is used tocarry out simulations of chemical processes with stan-dard equipments, gPROMS allows custom mathemati-cal description of new systems by writing modelequations within the simulation software itself. Oncethe new equipment model has been implemented, itcan be eventually connected with other units (e.g.pumps, measuring instruments, etc.), performingsimulation of the entire plant.

Several useful tools are also contained in gPROMSto help the user during model building: in particular,all information describing the system (variables,parameters and equations) must be declared in differ-ent entities following a hierarchical sequence (asshown in Fig. 5). In this way, a highly structuredmodel is built, which allows the analysis of systembehaviour under different conditions.

In the present case, the model has been imple-mented in two hierarchical levels (cell pair and stack).Input data used in simulation are summarised inTable 1.

Finally, the x-domain has been discretised in 50elements, each one with a length dx= 2mm, thus pro-viding information on variables along the wholecompartments.

3. Experimental

3.1. Experimental set-up and procedure

In order to validate the proposed model, experi-mental measurements of voltage and power densitywere carried out on a lab-scale stack, at FUJIFILM lab-oratories in Tilburg (the Netherlands), using differenttype of commercial spacers. A simplified scheme ofthe experimental apparatus is shown in Fig. 6. Theadopted equipment has a standard design for electro-dialysis process (Deukum GmbH, Germany), with

electrodes of DSA-type. Membranes are FUJIFILMmanufactured ion-exchange membranes with 10cm� 10 cm nominal dimensions (open area); other rel-evant physical properties of IEMs are listed in Table 2.Woven spacers with two different nominal thickness(365 and 450 lm) were used in order to investigate theeffect of channel thickness on process efficiency andto validate the model at different conditions. Twoperistaltic pumps (Hosepump Masterflex PW, Burt

Fig. 5. Model structure within gPROMS� Model Builder.

Table 1Summary of input data implemented in gPROMS

Input variables

Seawater inlet concentrations

Brine inlet concentrations

Seawater feed flow rate

Brine feed flow rate

Parameters

Membrane properties

CEM permselectivity

AEM permselectivity

CEM areal resistance

AEM areal resistance

CEM/AEM thickness

Cell geometry

Cell pair length and width

Seawater/brine compartments thickness


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

Process Equipment Inc., USA) were used to allowcirculation of both solutions through the stack; anadditional pump of the same type was used for circu-lation of electrode rinse solution.

Measurements were carried out using 0.5M NaClaqueous solution as “artificial seawater” and 5.4MNaCl aqueous solution as “artificial brine”. Anaqueous solution of K3Fe(CN)6/K4Fe(CN)6 0.1M andNaCl 2.5M was used as electrode rinse solution: theaddition of NaCl (with a concentration between 0.5and 5.4M) allowed the reduction of concentration gra-dient between electrode/central compartments, aswell as the improvement of conductivity of the elec-trode rinse solution itself.

The stack was connected to feed hydraulic circuitsby manifolds and to the measuring instrument (Auto-lab potentiostat/galvanostat PGSTAT100, Metrohm,USA) through electrical terminals connected with elec-trodes. A picture of experimental set-up is shown inFig. 7, in which it is possible to see the cells stack, aswell as the other relevant elements of the system.

Experimental measurements consisted of standardcyclic voltammetric analyses in galvanostatic mode[16]: in particular, each measurement was carried outimposing a ramp for the current from 0 to �0.2A with

a step of �0.5mA/s; the measuring instrumentrecords the potential difference (stack voltage) at ter-minals (Fig. 8). Each proof was repeated three timesin order to guarantee reproducibility of the measure-ment. Once the stack voltage is measured, the stack

Fig. 7. Experimental set-up.(1) Potentiostat/galvanostat; (2) cells stack; (3) pumps; (4)electrode rinse solution tank and (5) seawater/brineinlet–outlet tanks.

Brine IN

CELLSSTACK

Potentiostat/Galvanostat

Brine INtank

Seawater INtank

Seawater OUTtank

Brine OUTtank

Electrode RinseSolution

tank

Seawater IN

ERS IN ERS OUT

Seawater OUT

Brine OUT

ERS

Fig. 6. Scheme of experimental set-up.

Table 2Physical properties of used membranes

Membrane Permselectivitya (%) Areal resistance (X cm2) Thickness (lm)

CEM 90 2.6 120

AEM 65 1.1 120

aValues measured with NaCl solutions 0.5M for diluate and 4M for concentrate.


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

resistance can be easily calculated regarding Ohm’slaw (Eq. (19)). In fact, Rstack represents the slope of theexperimental curve in the graph Estack/I presented inFig. 8, while the intercept with the y-axis representsOCV for the system. Finally, electric power was

calculated as the product of stack voltage and appliedcurrent, and power density was estimated normalisingthe latter with respect to membrane area and cellsnumber.

0.00

0.10

0.20

0.30

0.40

0.50

-0.20 -0.15 -0.10 -0.05 0.00

Pote

ntia

l [V]

Electric current [A]

4 cells

OCV

Fig. 8. Ciclic voltammetry galvanostatic on a stackequipped with four cells.Stack equipped with 450lm spacers in electrode/centralcompartments; flow rate in single channel: 13mL/min.

0.0

0.5

1.0

1.5

2.0

0 2 4 6 8 10 12 14

Res

ista

nce

[Ω]

No. cells

Exp measurements

Blank measurement

Rblank

Fig. 9. Experimental measurements of stack and blankresistance.Exp measurements with 4-/6-/9-/12-cells stack equippedwith 450lm spacers in electrode/central compartments;flow rate in single channel: 13mL/min. Blankmeasurement is related to a stack equipped only with 1CEM and 2 end-spacers.

0.00

0.05

0.10

0.15

0.20

0.0 0.2 0.4 0.6 0.8 1.0

Pow

er O

utpu

t, P

[W]

Stack Voltage, Estack [V]

12

9

6

4

No. of cells

Fig. 10. Electric power for SGP-RE stack at differentnumber of cells.Exp measurements with 4-/6-/9-/12-cells stack equippedwith 450 lm spacers in electrode/central compartments;flow rate in single channel: 13mL/min.


0.00

0.05

0.10

0.15

0.20

0.0 0.2 0.4 0.6 0.8 1.0

"Cor

rect

ed"

Pow

er ,

P corr [W

]

12

9

6

4

No. of cells

Fig. 11. “Corrected” power for SGP-RE stack at differentnumber of cells.Exp measurements with 4-/6-/9-/12-cells stack equippedwith 450 lm spacers in electrode/central compartments;flow rate in single channel: 13mL/min.


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

3.2. Experimental measurements

3.2.1. Electrode compartments (blank) resistance

Experimental measurements of stack voltage areclearly affected also by the resistance of electrodecompartment (Rblank). In order to be able to correctthe measurements for obtaining the generic cell pairvoltage (i.e. the voltage theoretically obtainable with anegligible blank resistance), Rblank must be estimated.Evaluation of blank resistance was performed fromthe experimental measurements of stack resistancewith different number of cell pairs (Fig. 9). In fact, asindicated by Eq. (18), Rblank can be considered as theresistance of a “0-cells” stack and therefore can beevaluated by extrapolation of experimental data in agraph Rstack/N.

Measurements were carried out using stacksequipped with 450lm spacers and with 4–6–9–12 cellpairs. Feed solutions flow rate was regulated in orderto maintain a constant flow rate in single channelbetween 11 and 13mL/min (corresponding to a linearvelocity close to 1 cm/s). Results are shown in Fig. 9.

Blank resistance was also directly measured carry-ing out a cyclic voltammetry on a stack equipped onlywith the two end-spacers and a single CEM betweenthem. The latter is added to avoid the passage of elec-trode rinse solution directly from one compartment tothe other, allowing the flux of solution only throughexternal manifolds. Therefore in this case, the systemdoes not contain central compartments, but just theelectrode ones. Result of this experimental measure-ment is also reported in Fig. 9 for comparison.

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0 0.1 0.2 0.3 0.4

Pow

er D

ensi

ty, P

d [W

/m2 ]

Pow

er D

ensi

ty, P

d [W

/m2 ]

Pow

er D

ensi

ty, P

d [W

/m2 ]

Pow

er D

ensi

ty, P

d [W

/m2 ]

4 cells

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0 0.2 0.4 0.6

6 cells

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0 0.2 0.4 0.6 0.8

Stack Voltage, Estack [V] Stack Voltage, Estack [V]

Stack Voltage, Estack [V] Stack Voltage, Estack [V]

9 cells

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0 0.2 0.4 0.6 0.8 1.0 1.2

12 cells

Fig. 12. Model validation for cells stack equipped with 450lm spacers at different number of cell pairs.N Experimental data; s “Corrected” data; Model predictions: – · –·· b= 0.7; ____ b= 0.8; Model predictions for 200 cellsstack: ........ b= 0.7; _ _ _ b= 0.8.


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

Fig. 9 shows a clearly linear dependence of stackresistance on cells number, thus confirming theoreticalexpectations. On the other side, blank resistance eval-uated as intercept of the regression line with they-axis results roughly 10% higher than the experimen-tal value directly measured. Such discrepancy can beattributed to the different configuration of the stack inthe blank measurement: in fact, when the stack isequipped with a generic number of cells, the passageof solution between electrode/central compartmentscannot be completely avoided; while, in the case ofblank measurement, the stack has only a single CEMbathed to both electrode compartments and sucheffect is absent. Because of this low discrepancy, theextrapolated value from regression line was chosen asmeasure of Rblank for subsequent calculations andmodel validation.

3.2.2. Electric power and power density

The power output and efficiency for an “ideal”single cell pair, i.e. not affected by the resistance ofelectrode compartments, has been evaluated by a“corrected” stack voltage, calculated as the product ofthe overall cells resistance and the current. The cellsresistance was calculated from the stack resistance(Fig. 9) subtracting the contribution of the blank. Inthis way, it is possible to calculate an electric poweroutput, which is not affected by the resistance of theelectrode compartments, thus more consistent withthe definition given in Eq. (21). Relevant results for 4–6–9–12 cells stack are shown in Figs. 10 and 11. It isworth mentioning that experimental range of currentdensity was rather limited, thus allowing the construc-tion of a small portion of the quadratic trend reportedin the figure. On the other side, given the importance

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0 0.2 0.4 0.6 0.8

Pow

er D

ensi

ty, P

d [W

/m2 ]


9 cells

Fig. 13. Model validation for cells stack equipped with365lm spacers.NExperimental data; s “Corrected” data; Model pre-dictions for 9 cells stack: –· –·· b= 0.7; ____ b= 0.8; Modelpredictions for 200 cells stack: ........ b= 0.7; _ _ _ b= 0.8.

0.00

0.05

0.10

0.15

0.20

0 10 20 30 40 50

flow

rate

in c

hann

el, Q

i [x1

0-6m

3 /s]

cell pair position in the stack

s = 0.2 mms = 0.5 mms = 1.0 mm

Fig. 14. Estimated flowrates distribution for a 50-cellsstack.CFD simulation of a 50-cells stack with SISO configurationat different distributor thickness. Total feed flow rate:Qtot = 4� 10�6m3/s; spacer thickness of seawater/brinecompartments: d= 200 lm [17].

Table 3Effect of non-uniform flow rates distribution. Summary of case studies. Spacer thickness of seawater/brinecompartments: d=200 lm [17]

s1 = 0.2mm s2 = 0.5mm s3 = 1.0mm

Qtot,A = 4� 10�6m3/s (240mL/min) CASE # A1 CASE # A2 CASE # A3

Qtot,B = 2� 10�5m3/s (1,200mL/min) CASE # B1 CASE # B2 CASE # B3

Qtot,C = 4� 10�5m3/s (2,400mL/min) CASE # C1 CASE # C2 CASE # C3


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

of plotting the entire curve in order to evaluate themaximum power output of each configuration, in Sec-tion 4 the entire curve predicted by the model is plot-ted although comparison with experiments is limitedto the presented data.

As expected, power output increases whenincreasing the number of cells (Fig. 10). Moreover,Fig. 11 shows how the “corrected” values of powerdensity are significantly higher than measuredvalues, thus indicating a strong influence of the

0.00

0.02

0.04

0.06

0.08

0.10

0.00 0.02 0.04 0.06 0.08 0.10

Cel

l Pai

r Vol

tage

, Ece

ll [V]

Cel

l Pai

r Vol

tage

, Ece

ll [V]

Cel

l Pai

r Vol

tage

, Ece

ll [V]

flow direction, x [m]

cell pair 1cell pair 10cell pair 15cell pair 20

case A1

0.0

0.5

1.0

1.5

2.0

0.00

0.02

0.04

0.06

0.08

0.10

0 10 20 30 40 50cell pair position in the stack

Qi

Estack

case A1

0.00

0.02

0.04

0.06

0.08

0.10

0.00 0.02 0.04 0.06 0.08 0.10



case A2

0.0

0.5

1.0

1.5

2.0

0.00

0.02

0.04

0.06

0.08

0.10

0 10 20 30 40 50


case A2

Qi

Estack

0.00

0.02

0.04

0.06

0.08

0.10

0.00 0.02 0.04 0.06 0.08 0.10



case A3

0.0

0.5

1.0

1.5

2.0

0.00

0.02

0.04

0.06

0.08

0.10

0 10 20 30 40 50

flow

rate

in c

hann

el, Q

i [x1

0-6 m

3 /s]

flow

rate

in c

hann

el, Q

i [x1

0-6 m

3 /s]

flow

rate

in c

hann

el, Q

i [x1

0-6 m

3 /s]

aver

age

Stac

k Vo

ltage

, Est

ack [

V]av

erag

e St

ack

Volta

ge, E

stac

k [V]

aver

age

Stac

k Vo

ltage

, Est

ack [

V]


case A3

Qi

Estack

Fig. 15. Effect of non-uniform flowrates distribution on stack voltage (cases A).


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

blank resistance on the power output especially forstacks with low number of cells.

4. Model results and validation

4.1. Validation procedure and results

A first model tuning/validation has beenperformed by comparing model predictions andexperimental measurements, fixing the operating con-ditions adopted during experiments and alreadydescribed in the previous section. In particular, inorder to reproduce the tests with variable current den-sity, simulations were performed varying the stackinternal current density by acting on the external load(Ru), thus reproducing the same current ramp experi-mentally imposed by the galvanostat.

Moreover, in order to neglect the effect of blankresistance on total stack resistance and to compare the“corrected power density” values experimentallyworked out, simulations of 100, 150 and 200 cell pairsstacks were carried out, finding that already with 200cell pairs bank resistance was negligible.

A tuning of the adjustable parameter b was alsoperformed, and simulations with values of b ¼ 0:7and b ¼ 0:8 will be shown in the figures.

Model validation results are shown in Fig. 12 forcells stacks equipped with FUJIFILM membranes(specs given in Table 2) and 450 lm woven spacers;the four different cases with 4–6–9–12 cells within the

stack have been considered, where flow rates are uni-formly distributed in each cell.

Model predictions follow fairly well experimentaltrends, also indicating the parabolic shape of thepower output vs. stack voltage theoretically foreseen.

As physically expected, estimated stack voltage ishigher when permselectivity correction factor isincreased. However, for the investigated cases, thevalue of b between 0.70 and 0.80 seems to be able topredict fairly well the system behaviour, with theexception of the case of 12 cell pairs, in which theoptimum value for the model parameter seems to becloser to 0.7. However, this discrepancy could be dueto non-ideal behaviour of the system, related forexample to the presence of non-uniform flow ratesdistribution and, more importantly, short-cut currentswithin the stack.

It is worthnoting how such values of b are quitelow (the expected reduction of permselectivity passingfrom 4 to 5M solution is in the range of 5–10%), thusindicating that such correction factor likely accountfor other non-ideality factors such as the shadoweffect of spacer on membrane, which can undoubtedlyreduce the effective membrane area for ion passage.

Validation of the model was carried out also for365 lm spacers (Fig. 13), and similar considerationsapply to this case.

4.2. Effect of non-uniform flow rates distribution: modelcoupling with CFD simulations [17]

When the assumption of ideal flow rates distribu-tion through the stack is not acceptable anymore, e.g.in the case of a very thin distribution channel, theneed for assessing the effect of the real non-uniformdistribution arises and requires a more indepth math-ematical description of the problem also taking intoaccount the fluid flow behaviour of the system, bothin the distribution channels and in each single com-partment. This is part of a detailed analysis, whichhas been performed using computational fluiddynamics (CFD) tools, widely discussed in anotherwork [17].

In this section, the effect of non-uniform flowratesdistribution on process efficiency is discussed. Resultsof CFD simulations of a 50-cells stack with a single-input single-output (SISO) configuration [17] areshown in Fig. 14, showing the flow rates distributioninside the stack using different thicknesses for the dis-tribution channel (s).

CFD results were implemented within the presentmodel by adding an equation relating single channelflow rate to the “normalised number of cell pair”,

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.00 0.02 0.04 0.06 0.08 0.10

Pow

er D

ensi

ty, P

d [W

/m2 ]


case A1case A2case A3

Pd,ave= 1.96W/m 2

Pd,ave= 1.99 W/m 2

Fig. 16. Effect of non-uniform flowrates distribution onpower density (cases A).


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

defined as the id. number of the cell pair divided bythe total number of cells (N). Such equation has beenfound to have a quadratic form [17] and relevant coef-

ficients were related to the geometrical and operatingconditions examined, namely the thickness of distribu-tion channel, s, and total feed flow rate, Qtot (Table 3).

0.00

0.02

0.04

0.06

0.08

0.10

0.00 0.02 0.04 0.06 0.08 0.10

Cel

l Pai

r Vol

tage

, Ece

ll [V]

Cel

l Pai

r Vol

tage

, Ece

ll [V]

Cel

l Pai

r Vol

tage

, Ece

ll [V]



case B1

0.0

0.5

1.0

1.5

2.0

0.00

0.02

0.04

0.06

0.08

0.10

0 10 20 30 40 50

flow

rate

in c

hann

el, Q

i [x1

0-6 m

3 /s]

flow

rate

in c

hann

el, Q

i [x1

0-6 m

3 /s]

flow

rate

in c

hann

el, Q

i [x1

0-6 m

3 /s]

aver

age

Stac

k Vo

ltage

, Est

ack [

V]av

erag

e St

ack

Volta

ge, E

stac

k [V]

aver

age

Stac

k Vo

ltage

, Est

ack [

V]


case B1

Qi

Estack

0.00

0.02

0.04

0.06

0.08

0.10

0.00 0.02 0.04 0.06 0.08 0.10



case B2

0.0

0.5

1.0

1.5

2.0

0.00

0.02

0.04

0.06

0.08

0.10

0 10 20 30 40 50


case B2

Qi

Estack

0.00

0.02

0.04

0.06

0.08

0.10

0.00 0.02 0.04 0.06 0.08 0.10



case B3

0.0

0.5

1.0

1.5

2.0

0.00

0.02

0.04

0.06

0.08

0.10

0 10 20 30 40 50


case B3

Qi

Estack

Fig. 17. Effect of non-uniform flowrates distribution on stack voltage (cases B).


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

Nine case studies were analysed in order to inves-tigate the effect of different flowrates distribution,generated by the three possible values for distributorthickness and total feed flow rate, on process operat-ing and performance parameters (Table 3).

Simulations results related to cases “A”(Qtot,A = 4� 10�6m3/s) are shown in Fig. 15. Note thatin case A1, where the distributor thickness is thinner,the average stack voltage is clearly minimised for thecell pairs in the middle of the stack: this is a conse-quence of the lower flowrate within these cells, i.e. ofthe higher residence time within compartments. Thiseffect is rather negligible in the other cases (A2, A3),where a thicker distributor allows more uniform flowrates along the stack [17].

The effect of flow rates distribution on obtainablepower density is shown in Fig. 16. Comparing theestimated value of power density for the cases A1–2–3, it can be seen in case A1 (thinner distributor andnon-uniform flow rates distribution) a reduction byroughly 5% of the mean power density with respectto cases A2 and A3. At the contrary, for cases A2 andA3 the obtainable power density is nearly the same,due to the very similar flow rates distributions.

A very similar outcome can be observed inFigs. 17–20, were results for cases B and C are shown.Also in this case, the configurations with the thinnerdistributor channel present a minimum in stack poten-tial for the central cell pairs and, consequently, areduction in the output power density, which is, how-ever, negligible in all cases. This is likely related alsoto the very high flow rates in the stack, which keep

the driving force in all cell pairs sufficiently high notto affect significantly the overall process performance.

4.2.1. Estimating the obtainable net power

Another important goal achievable through themerging of CFD and process modelling approach hasbeen the estimation of the obtainable electric powerfrom reverse electro-dialysis (RED) system, once therequired pumping power has been taking into account.In order to do this, net power is calculated for all ofinvestigated case studies (Table 4), using Eqs. (22) and(23). Results are listed in Table 4. Pressure drop inTable 4 are calculated from CFD simulation and dis-cussed in another work [17] and, although they mayrefer to an ideal stack geometry, their use can shownof the multiscale approach can be usefully adopted forthe simulation of overall process performance.

In particular, Table 4 shows how increasing flowrates generally increases the output power density,but at the same time it can dramatically increasepumping power, thus leading to a negative influenceon the net power density. As an example, in case C1,(higher total feed flow rate and thinner distributorthickness), the pressure drop within the system is sohigh that pumping power is eventually higher thanobtainable electric power, thus indicating that an opti-mal compromise has to be chosen for the overall pro-cess optimisation in terms of both geometrical andoperating conditions.

4.3. Influence of membrane properties on power density:perspectives for achieving the goal of 8W/m2 powergeneration

A “must” condition for the technological break-through of the SGP-RE process is the achievementof a power density sufficiently high to allow for theproduction of electricity at a cost competitive withother renewable energy technologies. With thisregard, one of the main goals of the REAPower pro-ject is the development of very thin IEMs in orderto reduce dramatically the stack electrical resistance,which is already significantly reduced by the use ofhighly conductive solutions (seawater instead ofriver water is used as dilute solution). Increasingand keeping high values of permselectivity for thinmembranes is also another important goal, whichwould allow a remarkable enhancement of the SGP-RE process performance.

On this basis, the model has been adopted forpredicting the achievable power density assuming adecrease in membrane thickness (and, proportion-

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.00 0.02 0.04 0.06 0.08 0.10

Pow

er D

ensi

ty, P

d [W

/m2 ]


case B1case B2case B3

Pd,ave= 2.36 W/m2

Pd,ave= 2.37 W/m2

Fig. 18. Effect of non-uniform flowrates distribution onpower density (cases B).


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

ally, also of the IEMs resistance) from the previous120lm to values of 20 lm and fixing the AEMpermselectivity to the values of 0.65 (as in the

standard case shown in previous simulations) and0.85, which is the target value expected to bereached within the project. The same geometry

0.00

0.02

0.04

0.06

0.08

0.10

0.00 0.02 0.04 0.06 0.08 0.10

Cel

l Pai

r Vol

tage

, Ece

ll [V]

Cel

l Pai

r Vol

tage

, Ece

ll [V]

Cel

l Pai

r Vol

tage

, Ece

ll [V]



case C1

0.0

0.5

1.0

1.5

2.0

0.00

0.02

0.04

0.06

0.08

0.10

0 10 20 30 40 50

flow

rate

in c

hann

el, Q

i [x1

0-6 m

3 /s]

flow

rate

in c

hann

el, Q

i [x1

0-6 m

3 /s]

flow

rate

in c

hann

el, Q

i [x1

0-6 m

3 /s]


case C1

Qi

Estack

0.00

0.02

0.04

0.06

0.08

0.10

0.00 0.02 0.04 0.06 0.08 0.10



case C2

0.0

0.5

1.0

1.5

2.0

0.00

0.02

0.04

0.06

0.08

0.10

0 10 20 30 40 50


case C2

Qi

Estack

0.00

0.02

0.04

0.06

0.08

0.10

0.00 0.02 0.04 0.06 0.08 0.10



case C3

0.0

0.5

1.0

1.5

2.0

0.00

0.02

0.04

0.06

0.08

0.10

0 10 20 30 40 50

aver

age

Stac

k Vo

ltage

, Est

ack [

V]av

erag

e St

ack

Volta

ge, E

stac

k [V]

aver

age

Stac

k Vo

ltage

, Est

ack [

V]


case C3

Qi

Estack

Fig. 19. Effect of non-uniform flowrates distribution on stack voltage (cases C).


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

adopted in the previous paragraph was used, i.e.spacer thickness of seawater/brine compartmentsd= 200lm and cell pair dimensions 10� 10 cm. Inorder to take into accountof the effect of high salin-ity and spacer shadow on membrane surface, bvalue is kept at 0.75, as already done during themodel tuning/validation procedure, thus guarantee-ing a conservative estimation of predicted values.

A first important remark on the effect of mem-brane thickness on the stack resistance is IEMs resis-tance that represents a significant portion of theoverall resistance when membrane thickness is around100lm, while for thinner membranes seawater

compartment resistance may play a dominant role incontrolling the overall cell pair resistance. In thisrespect, a crucial aspect will be the choice of compart-ment thicknesses. In fact, a thicker brine compartmentwould allow minor problems of fouling and pluggingwithout significantly affecting the stack resistance, onthe other side, seawater compartment will need athinner spacer, thus requiring a better control of foul-ing and plugging by a suitable pretreatment of feedsolution (see Fig. 21).

Finally, power density predictions are shown inFig. 22, which indicates how reducing IEMs thicknesscan significantly enhance the power density, up tovalues of about 6.7W/m2 for the 20lm IEMs and 0.65AEM permselectivity. Moreover, increasing the AEMpermselectivity to the target value of 0.85 would allowfor a further raise in power density with a maximumexpected value of more than 8.5W/m2, thus indicat-ing a tremendous potential for improvement of thepresented technology.

5. Conclusions

SGP-RE has been recently recognised as a promis-ing technology for energy generation from salinitygradients. Latest interests have also arisen with regardto the exploitation of salinity gradients from brinesand seawater, which can guarantee much higher driv-ing forces to the SGP process.

Aim of the present work, carried out within theEU-FP7-funded REAPower project, has been thedevelopment of a mathematical model for the SGP-REprocess using seawater and concentrated brine.

Starting from a literature review, the model wasdeveloped focusing on its capability to simulate the

Table 4Estimated net power for case studies of Table 3

Distributorthickness

Feed flow rate Pressuredrop

Pumpingpower

Powerdensity

Poweroutput

Netpower

Powerlosses

CASE#

s(mm)

Qtot

(�10�6m3/s)DP(Pa)

Ppump

(W)Pd

(W/m2)Pout

(W)Pnet

(W)Ploss

(%)

A1 0.2 4 1,326 0.015 1.96 0.98 0.97 1.5

A2 0.5 4 391 0.004 1.99 1.00 0.99 0.4

A3 1.0 4 330 0.004 1.99 1.00 0.99 0.4

B1 0.2 20 7,118 0.407 2.36 1.18 0.77 34.4

B2 0.5 20 2,054 0.117 2.37 1.19 1.07 9.9

B3 1.0 20 1,706 0.097 2.37 1.19 1.09 8.2

C1 0.2 40 15,736 1.798 2.42 1.21 �0.59 148.7

C2 0.5 40 4,410 0.504 2.43 1.21 0.71 41.5

C3 1.0 40 3,567 0.408 2.43 1.22 0.81 33.5

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.00 0.02 0.04 0.06 0.08 0.10

Pow

er D

ensi

ty, P

d [W

/m2 ]


case C1case C2case C3

Pd,ave= 2.42 W/m 2

Pd,ave= 2.43 W/m2

Fig. 20. Effect of non-uniform flowrates distribution onpower density (cases C).


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

system behaviour when concentrated salt solutionswere adopted. A multiscale approach has been used,leading to a multihierarchical model able to predict

the effects of operating conditions and geometricalfeatures on the process performance, even allowingthe simulation of non-uniform flow rates distributionwithin a 50-cells stack to the coupling with CFDsimulations. The developed model was tuned andvalidated by comparison with experimental informa-tion purposely collected on a laboratory test rig.

Finally, some technology perspectives were alsogiven by performing simulations aiming at the predic-tion of power density achievable when using thin andselective ionic exchange membranes, indicating thefeasibility of reaching a target power density of morethan 8.5W/m2 of cell pair.

Nomenclature

A –– membrane area (m2)

b –– membrane width (m)

B´ –– parameter of Islam et al.’ equation(m1/2mol�1/2)

B01 –– parameter of Islam et al.’ equation

(Sm3mol�3/2)

B02 –– parameter of Islam et al.’ equation

(m3/2mol�1/2)

Bc –– second virial coefficient of Pitzerequation (kgmol�1)

c –– salt concentration (molm�3)

Cc –– third virial coefficient of Pitzerequation (Kg2mol�2)

DNaCl –– co-ion diffusion coefficient (m2 s�1)

Dw –– water diffusion coefficient (m2 s�1)

Estack –– stack voltage (V)

F –– Faraday constant (96,490Cmol�1)

f –– obstruction factor (–)

F´ –– parameters of Islam et al.’ equation (–)

fc –– first virial coefficient of Pitzer equation(–)

I –– electric current (A)

j –– current density (Am�2)

J0w –– volumetric water flux (m3m�2 s�1)

Jcoul, Jcit –– counter-ions/co-ions molar flux(molm�2 s�1)

Jtot –– salt molar flux (molm�2 s�1)

Jw –– water flux (molm�2 s�1)

L –– channel length (m)

m –– electrolyte molal concentration(mol kg�1)

MH2O –– molar weight of water(18� 10�3 kgmol�1)

N –– number or cell pairs (–)

Pd –– power density (Wm�2)

Ploss –– power loss (%)

Pnet –– net power output (W)

52.0% 39.6% 32.0% 26.9% 20.4%

13.0% 9.5% 7.6% 6.3% 4.8%24.6%

35.7%42.4%

46.9%

52.6%

10.4%

15.1%

17.9%

19.8%

22.3%

0.0

1.0

2.0

3.0

4.0

5.0

6.0

20 40 60 80 120

Are

al R

esis

tanc

e [x

10-4 Ω

·m2 ]

AEM/CEM thickness [μm]

AEM

CEM

brine

seawater

Fig. 21. Influence of IEMs thickness on resistance of thesystem. Simulation of a 1,000 cells stack assuming a lineardecreasing in IEMs resistance with IEMs thickness.aAEM= 0.65, aCEM= 0.90. Spacer thickness of seawater/brine compartments d= 200lm.

0.0

2.0

4.0

6.0

8.0

10.0

0 20 40 60 80 100 120 140

Mea

n Po

wer

Den

sity

[W/m

2 ]

AEM/CEM thickness [μm]

αAEM

Fig. 22. Effect of IEMs thickness on power density.Simulation of a 1,000 cells stack assuming a lineardecreasing in IEMs resistance with IEMs thickness.N: aAEM= 0.65, aCEM= 0.90 ; s: aAEM= 0.85, aCEM=0.90.


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

Pout –– power output (W)

Ppump –– pumping power (W)

Qb, Qs –– brine/seawater flow rate in singlechannel (m3 s�1)

Qtot –– total feed flow rate (m3 s�1)

R –– universal gas constant (8.314 Jmol�1

K�1)

RAEM, RCEM –– AEM/CEM areal resistance (Xm2)

Rb, Rs –– brine/seawater areal resistance (Xm2)

Rblank –– electrode compartments (blank)resistance (Xm2)

Rcell –– cell pair resistance (Xm2)

Rstack –– cells stack resistance (Xm2)

Ru –– external load (Xm2)

T –– temperature (K)

Ta, Tc –– transport number of cation and anionin solution (–)

TCEMc –– transport number of cation into CEM

(–)

x –– flow direction (m)

z –– ion valence (–)

Greek letters

aAEM, aCEM –– AEM/CEM permselectivity (–)

b –– permselectivity correction factor (–)

c± –– mean activity coefficient (–)

db, ds –– brine/seawater compartment thickness(m)

dm –– membrane thickness (m)

DP –– pressure drop (Pa)

e –– dielectric constant (–)

g –– viscosity of solution (Pa s)

gp –– pump efficiency (–)

K –– equivalent conductivity (Sm2mol�1)

K0 –– equivalent conductivity at infinitedilution (126.5� 10�4 Sm2mol�1)

qH2O –– water density (1,000 kgm�3)

Acknowledgements

This work has been performed within the REA-Power (Reverse Electro dialysis Alternative Power

production) project [2], funded by the EU-FP7 pro-gramme (Project No. 256736). The authors are verygrateful to Winod Bhikhi and Jacko Hessing of FUJI-FILM Tilburg Research Laboratories TRL for theirsupport within experimental activities.

References

[1] R.E. Pattle, Production of electric power by mixing fresh andsalt water in the hydroelectric pile, Nature 174 (1954) 660.

[2] REAPower, http://reapower.eu/.[3] H. Strathmann, Ion-Exchange Membrane Separation Pro-

cesses, Membrane Science and Technology, vol. 9, Elsevier,Amsterdam, 2004, pp. 1–348.

[4] J.N. Weinstein, F.B. Leitz, Electric power from differences insalinity: The dialytic battery, Science 191 (1976) 557–559.

[5] C. Forgacs, R.N. O’Brien, Utilization of membrane processesin the development of non-conventional renewable energysources, Chem. Can. 31 (1979) 19–21.

[6] R.E. Lacey, Energy by reverse electrodialysis, Ocean Eng. 7(1980) 1–47.

[7] E. Brauns, Salinity gradient power by reverse electrodialysis:effect of model parameters on electrical power output, Desali-nation 237 (2009) 379–391.

[8] J. Veerman, J.W. Post, M. Saakes, S.J. Metz, G.J. Harmsen,Reverse electrodialysis: A validated process model for designand optimization, Chem. Eng. J. 166 (2011) 256–268.

[9] J.O.M. Bockris, A.K.N. Reddy, Modern Electrochemistry,vol. 1 (Ionics), Plenum Press, New York, NY, 1998.

[10] K.S. Pitzer, Thermodynamics of electrolytes. I. theoreticalbasis and general equations, J. Phys. Chem. 77 (1973)268–277.

[11] National Physical Laboratory UK, http://www.kayelaby.npl.co.uk/chemistry/3_9/3_9_6.html (accessed January, 2012).

[12] S.S. Islam, R.L. Gupta, K. Ismail, Extension of the falkenha-gen-leistkelbg to the electrical conductance of concentratedaqueous electrolytes, J. Chem. Eng. Data 36 (1991) 102–104.

[13] A. de Diego, A. Usobiaga, J.M. Madariaga, Critical compari-son among equations derived from the falkenhagen model tofit conconductimetric data of concentrated electrolyte solu-tions, J. Electroanal. Chem. 446 (1998) 177–187.

[14] M. Turek, B. Bandura, Renewable energy by reverse electrodi-alysis, Desalination 205 (2007) 67–74.

[15] gPROMS Model Builder 3.3.1 (2010).[16] J. Veerman, R.M. de Jong, M. Saakes, S.J. Metz, G.J. Harmsen,

Reverse electrodialysis: Comparison of six commercial mem-brane pairs on the thermodynamic efficiency and power den-sity., J. Membr. Sci. 343 (2009) 7–15.

[17] L. Gurreri, A. Tamburin, A. Cipollina, G. Micale, CFD analysisof the fluid flow behavior in a reverse electrodialysis stack,Desalination and Water Treatment, 2012, in press. doi:10.1080/19443994.2012.705966.


Dow

nloa

ded

by [

Uni

vers

ita d

i Pal

erm

o], [

And

rea

Cip

ollin

a] a

t 10:

00 2

7 N

ovem

ber

2012

http://reapower.eu/

http://www.kayelaby.npl.co.uk/chemistry/3_9/3_9_6.html

http://www.kayelaby.npl.co.uk/chemistry/3_9/3_9_6.html

Date post:	23-Nov-2023
Category:	Documents
Upload:	independent
View:	0 times
Download:	0 times

Modelling the Reverse ElectroDialysis process with seawater and concentrated brines

Documents