Skip to content

Latest commit

 

History

History
368 lines (295 loc) · 26.2 KB

electrodialysis_1D.rst

File metadata and controls

368 lines (295 loc) · 26.2 KB
Raw

Electrodialysis (1D)

Introduction

Electrodialysis, an electrochemical separation technology, has been utilized to desalinate water for decades. Compared with other technologies, such as reverse osmosis (RO), electrodialysis shows advantages in desalinating brackish waters with moderate salinity by its less intense energy consumption, higher water recovery, and robust tolerance for adverse non-ionic components (e.g., silica and biological substances) in water. Between the electrodes of an electrodialysis stack (a reactor module), multiple anion- and cation-exchange membranes are alternately positioned and separated by spacers to form individual cells. When operated, electrodialysis converts electrical current to ion flux and, assisted by the opposite ion selectivity of cation- and anion-exchange membranes (cem and aem), moves ions from one cell to its adjacent cell in a cell-pair treatment unit (Figure 1). The ion-departing cell is called a diluate channel and the ion-entering cell a concentrate channel. Recovered (desalinated) water is collected from diluate channels of all cell pairs while the concentrate product can be disposed of as brine or retreated. More overview of the electrodialysis technology can be found in the References.

Figure 1. Schematic representation of an electrodialysis cell pair

Figure 1. Schematic representation of an electrodialysis cell pair

One cell pair in an electrodialysis stack can thus be treated as a modeling unit that can multiply to larger-scale systems. The presented electrodialysis 1D model establishes mathematical descriptions of ion and water transport in a cell pair and expands to simulate a stack with a specified cell-pair number. Modeled mass transfer mechanisms include electrical migration and diffusion of ions and osmosis and electroosmosis of water. By its name, "1D" suggests the variation of mass transfer and solution properties over length (l in Figure 1) is mathematically described. The model relies on the following key assumptions:

  • The concentrate and diluate channels have identical geometry.
  • For each channel, component fluxes in the bulk solution are uniform in b and s directions in Figure 1 but dependent on the l directions (the 1-dimensional assumption).
  • Steady state: all variables are independent of time.
  • Co-current flow operation.
  • Electrical current is operated below the limiting current.
  • Ideality assumptions: activity, osmotic, and van't Hoff coefficients are set at one.
  • All ion-exchange membrane properties (ion and water transport number, resistance, permeability) are constant.
  • Detailed concentration gradient effect at membrane-water interfaces is neglected.
  • Constant pressure and temperature through each channel.

Control Volumes

This model has two 1D control volumes for the concentrate and diluate channels.

  • Diluate
  • Concentrate

Ports

On the two control volumes, this model provides four ports (Pyomo notation in parenthesis):

  • inlet_diluate (inlet)
  • outlet_diluate (outlet)
  • inlet_concentrate (inlet)
  • outlet_concentrate (outlet)

Sets

This model can simulate electrodialysis for desalination of a water solution containing multiple species (neutral or ionic). All solution components ( H2O, neutral solutes, and ions) form a Pyomo set in the model. For a clear model demonstration, this document uses an NaCl water solution as an instance hereafter. The user can nevertheless expand the component set as needed to represent other feed water conditions.

This model can mathematically take a multi-component (i.e., > one salt molecule to be treated) as an input; nevertheless a multi-component solution creates unknown or difficult-to-specify parameters, e.g., the electrical transport numbers through membranes, the multi-ion diffusivity, etc., and physical relationships, which may result in ill-posed or ill-conditioned problems challenging the models' numerical solutions. While we continuously work on advancing our models to absorb new principles revealed by progressing research, we advise the users be very cautious with simulating multi-component systems by this programmed model for aspects stated above.

Table 1. List of Set
Description Symbol Indices
Time t [t] ([0])1
Length_domain x l × (0, 1) 2
Phase p ['Liq']
Component j ['H2O', 'Na+', 'Cl-']
Ion j ['Na+', 'Cl-'] 3
Membrane n/a ['cem', 'aem']
Notes

1 The time set index is set as [0] in this steady-state model and is reserved for the future extension to a dynamic model.

2 By the IDAES convention, the index of length_domain is normalized to a continuous set of (0, 1), which is discretized when differential equations in the model are solved by numerical methods such as "finite difference" discretization. In this documentation, x refers to the length dimension before normalization and carries a unit of [m].

3 "Ion" is a subset of "Component" and uses the same symbol j.

Degrees of Freedom

Applying this model to an NaCl solution yields 33 degrees of freedom (Table 2), among which temperature, pressure, and component molar flow rate are state variables that are fixed as initial conditions. The rest are parameters that should be provided in order to fully solve the model.

Table 2. List of Degree of Freedom (DOF)
Description Symbol Variable Name Index Units DOF Number 1
Temperature, inlet_diluate TD temperature None K 1
Temperature, inlet_concentrate TC temperature None K 1
Pressure, inlet_diluate pD temperature None Pa 1
Pressure, inlet_concentrate pC temperature None Pa 1
Component molar flow rate, inlet_diluate Nj, inD flow_mol_phase_comp [t], ['Liq'], ['H2O', 'Na+', 'Cl-'] mols − 1 3
Component molar flow rate, inlet_concentrate Nj, inC flow_mol_phase_comp [t], ['Liq'], ['H2O', 'Na+', 'Cl-'] mols − 1 3
Water transport number tw water_trans_number_membrane ['cem', 'aem'] dimensionless 2
Water permeability L water_permeability_membrane ['cem', 'aem'] m − 1s − 1Pa − 1 2
Voltage or Current 2 U or I voltage_applied or current_applied [t] V or A 1
Electrode areal resistance rel electrodes_resistance [t] Ωm2 1
Cell pair number n cell_pair_num None dimensionless 1
Current utilization coefficient ξ current_utilization None dimensionless 1
Channel height d channel_height none m 1
Membrane areal resistance r membrane_areal_resistance ['cem', 'aem'] Ωm2 2
Cell width b cell_width None m 1
Cell length l cell_length None m 1
Thickness of ion exchange membranes δ membrane_thickness ['cem', 'aem'] m 2
diffusivity of solute in the membrane phase D solute_diffusivity_membrane ['cem', 'aem'], ['Na+', 'Cl-'] m2s − 1 4
transport number of ions in the membrane phase tj ion_trans_number_membrane ['cem', 'aem'], ['Na+', 'Cl-'] dimensionless 4
Note

1 DOF number accounts for indices of corresponding variables.

2 A user should provide either current or voltage as the electrical input, in correspondence to the "Constant_Current" or "Constant_Voltage" treatment mode (configured in this model). The user also should provide an electrical magnitude that ensures an operational current below the limiting current of the feed solution.

Solution component information

To fully construct solution properties, users need to provide basic component information of the feed solution to use this model, including identity of all solute species (i.e., Na +, and Cl - for a NaCl solution), molecular weight of all component species (i.e., H2O, Na +, and Cl -), and charge and electrical mobility of all ionic species (i.e., Na +, and Cl -). This can be provided as a solution dictionary in the following format (instantiated by a NaCl solution).

ion_dict = {
         "solute_list": ["Na_+", "Cl_-"],
         "mw_data": {"H2O": 18e-3, "Na_+": 23e-3, "Cl_-": 35.5e-3},
         "electrical_mobility_data": {"Na_+": 5.19e-8, "Cl_-": 7.92e-8},
         "charge": {"Na_+": 1, "Cl_-": -1},
     }

This model, by default, uses H2O as the solvent of the feed solution.

Information regarding the property package this unit model relies on can be found here:

:pywatertap.property_models.ion_DSPMDE_prop_pack

Equations

This model solves mass balances of all solution components (H2O, Na +, and Cl - for an NaCl solution) on two control volumes (concentrate and diluate channels). Under the 1D treatment, balance equations are expressed as differential algebraic equations (DAE) when concerned variables are functions of length (x). The DAEs are solved in a discretization manner using the "finite difference" or "collocation" method implemented in Pyomo.DAE.

Mass balance equations are summarized in Table 3. Mass transfer mechanisms account for solute electrical migration, diffusion, water osmosis, and electroosmosis. Theoretical principles, e.g., continuity equation, Fick's law, and Ohm's law, to simulate these processes are well developed and some good summaries for the electrodialysis scenario can be found in the References.

Table 3 Mass Balance Equations
Description Equation Index set
Component mass balance $\left(\frac{\partial N_j (x)}{\partial x}\right)^{C\: or\: D}+J_j(x)^{C\: or\: D} b=0$ j ∈ [′H2O′,′Na+′,′Cl−′]
mass transfer flux, concentrate, solute $J_j^{C} = \left(t_j^{cem}-t_j^{aem} \right)\frac{\xi i(x)}{ z_j F}-\left(\frac{D_j^{cem}}{\delta ^{cem}} +\frac{D_j^{aem}}{\delta ^{aem}}\right)\left(c_j^C(x)-c_j^D(x) \right)$ j ∈ [′Na+′,′Cl−′]
mass transfer flux, diluate, solute $J_j^{D} = -\left(t_j^{cem}-t_j^{aem} \right)\frac{\xi i(x)}{ z_j F}+\left(\frac{D_j^{cem}}{\delta ^{cem}} +\frac{D_j^{aem}}{\delta ^{aem}}\right)\left(c_j^C(x)-c_j^D(x) \right)$ j ∈ [′Na+′,′Cl−′]
mass transfer flux, concentrate, H2O $J_j^{C} = \left(t_w^{cem}+t_w^{aem} \right)\frac{i(x)}{F}+\left(L^{cem}+L^{aem} \right)\left(p_{osm}^C(x)-p_{osm}^D(x) \right)\frac{\rho_w}{M_w}$ j ∈ [′H2O′]
mass transfer flux, diluate, H2O $J_j^{D} = -\left(t_w^{cem}+t_w^{aem} \right)\frac{i(x)}{F}-\left(L^{cem}+L^{aem} \right)\left(p_{osm}^C(x)-p_{osm}^D(x) \right)\frac{\rho_w}{M_w}$ j ∈ [′H2O′]

Additionally, several other equations are built to describe the electrochemical principles and electrodialysis performance.

Table 4 Electrical and Performance Equations
Description Equation
Electrical input condition $i(x) = \frac{I}{bl}$, for 'Constant_Current'; u(x) = U for 'Constant_Voltage'
Ohm's law u(x) = i(x)rtot(x)
Resistance calculation $r_{tot}(x)=n\left(r^{cem}+r^{aem}+\frac{d}{\kappa^C(x)}+\frac{d}{\kappa^D(x)}\right)+r_{el}$
Electrical power consumption P(x) = b∫0lu(x)i(x)dx
Water-production-specific power consumption $P_Q=\frac{P(x=l)}{3.6\times 10^6 nQ_{out}^D}$
Current efficiency for desalination $bi(x)\eta(x)=-\sum_{j \in[cation]}{\left[\left(\frac{\partial N_j ^D(x)}{\partial x}\right) z_j F\right]}$

All equations are coded as "constraints" (Pyomo). Isothermal and isobaric conditions apply.

Extended simulation

Membrane and interfacial potentials

This model supports extensive simulations of (1) the nonohmic potential across ion exchange membranes and (2) the Nernst diffusion layer. Users can customize these extenions via two configurations: has_nonohmic_potential_membrane that triggers the calculation of nonohmic potentials across ion exchange membranes and has_Nernst_diffusion_layer that triggers the simulation of a concentration-polarized Nernst diffusion layer including its ohmic and nonohmic potential changes. Based on a electrochemical cell setup in Figure 2 and established theoretical descriptions (References), our model accounts for the cross-membrane diffusion and Donnan potentials (nonohmic), ion concentration polarization in assumed Nernst diffusion layers (NDL), and the ohmic and nonohmic (i.e., diffusion) potentials across NDLs. These extensions make the model closer to the non-ideal physical conditions that can be encountered in real desalination practices.

Figure 2. Electrochemical cell setup for simulating Nernst diffusion layer and cross-membrane potential and concentration variations.

Figure 2. Electrochemical cell setup for simulating Nernst diffusion layer and cross-membrane potential and concentration variations.

Table 5 presents the equations underlying the two extensions assuming a 1:1 symmetric electrolyte such as NaCl.

Table 5 Essential equations supporting model extensions
Description Equation Condition
Limiting current density $i_{lim}(x) = i_{lim,0}\frac{c_b^D(x)}{c_b^D(0)}$ has_Nernst_diffusion_layer==True and limiting_current_density_method == LimitingCurrentDensityMethod.InitialValue
ilim(x) = AvBcbD(x) has_Nernst_diffusion_layer==True and limiting_current_density_method == LimitingCurrentDensityMethod.Empirical
$i_{lim}(x) = \frac{Sh F D_b c_b^D(x)}{d_H \left(t_+^{cem}-t_+ \right)}$ has_Nernst_diffusion_layer==True and limiting_current_density_method == LimitingCurrentDensityMethod.Theoretical
Nonohmic potential, membrane $\phi_m(x)=\frac{RT}{F} \left( t_+^{iem} - t_-^{iem} \right) \ln \left( \frac{c_s^R(x)}{c_s^L(x)} \right)$ has_nonohmic_potential_membrane == True
Ohmic potential, NDL $\phi_d^{ohm}(x)=\frac{FD_b}{\left(t_+^{iem}-t_+\right)\lambda}\ln\left(\frac{c_s^L(x)c_b^R(x)}{c_s^R(x)c_b^L(x)}\right)$ has_Nernst_diffusion_layer==True
Nonohmic potential, NDL $\phi_d^{nonohm}(x)=\frac{RT}{F}\left(t_+-t_-\right) \ln\left(\frac{c_s^L(x)c_b^R(x)}{c_s^R(x)c_b^L(x)}\right)$ has_Nernst_diffusion_layer==True
NDL thickness, cem $\Delta^{L/R}(x) = \frac{F D_b c_b^{L/R}(x)}{\left(t_+^{iem}-t_+ \right) i_{lim}(x)}$ has_Nernst_diffusion_layer==True
NDL thickness, aem $\Delta^{L/R}(x) = - \frac{F D_b c_b^{L/R}(x)}{\left(t_+^{iem}-t_+\right) i_{lim}(x)}$ has_Nernst_diffusion_layer==True
Concentration polarization ratio, cem $\frac{c_s^L(x)}{c_b^L(x)} = 1+\frac{i(x)}{i_{lim}(x)},\qquad \frac{c_s^R(x)}{c_b^R(x)} = 1-\frac{i(x)}{i_{lim}(x)}$ has_Nernst_diffusion_layer==True 1
Concentration polarization ratio, aem $\frac{c_s^L(x)}{c_b^L(x)} = 1-\frac{i(x)}{i_{lim}(x)},\qquad \frac{c_s^R(x)}{c_b^R(x)} = 1+\frac{i(x)}{i_{lim}(x)}$ has_Nernst_diffusion_layer==True

Note

1 When this configuration is turned off, ilim is considered as ∞ and the ratio becomes 1.

Some other modifications to previously defined equations are made to accommodate the two extensions. These are shown in Table 6.

Table 6 Other equation modifications under extensions
Original equation description Equation replacement Condition
Ohm's law u(x) = i(x)rtot(x) + ϕm(x) + ϕdohm(x) + ϕdnonohm(x) 1 has_nonohmic_potential_membrane == True and/or has_Nernst_diffusion_layer==True
Resistance calculation $r_{tot}(x)=n\left(r^{cem}+r^{aem}+\frac{d- \Delta_{cem}^L(x) - \Delta_{aem}^R(x)}{\kappa^C(x)}+\frac{d- \Delta_{cem}^R(x) - \Delta_{aem}^L(x)}{\kappa^D(x)}\right)+r_{el}$ has_Nernst_diffusion_layer==True
mass transfer flux, concentrate, solute $J_j^{C} = \left(t_j^{cem}-t_j^{aem} \right)\frac{\xi i(x)}{ z_j F}-\left(\frac{D_j^{cem}}{\delta ^{cem}}\left(c_{s,j}^{L,cem}(x)-c_{s,j}^{R,cem}(x) \right) +\frac{D_j^{aem}}{\delta ^{aem}} \left(c_{s,j}^{R,aem}(x)-c_{s,j}^{L,aem}(x) \right)\right)$ has_nonohmic_potential_membrane == True and/or has_Nernst_diffusion_layer==True
mass transfer flux, diluate, solute $J_j^{D} = -\left(t_j^{cem}-t_j^{aem} \right)\frac{\xi i(x)}{ z_j F}+\left(\frac{D_j^{cem}}{\delta ^{cem}}\left(c_{s,j}^{L,cem}(x)-c_{s,j}^{R,cem}(x) \right) +\frac{D_j^{aem}}{\delta ^{aem}} \left(c_{s,j}^{R,aem}(x)-c_{s,j}^{L,aem}(x) \right)\right)$ has_nonohmic_potential_membrane == True and/or has_Nernst_diffusion_layer==True
mass transfer flux, concentrate, H2O $J_j^{C} = \left(t_w^{cem}+t_w^{aem} \right)\frac{i(x)}{F}+\left(L^{cem} \left(p_{s, osm}^{cem, L}(x)-p_{s, osm}^{cem, R}(x) \right)+L^{aem} \left(p_{s, osm}^{aem, R}(x)-p_{s, osm}^{aem, L}(x) \right)\right)\frac{\rho_w}{M_w}$ has_Nernst_diffusion_layer==True
mass transfer flux, diluate, H2O $J_j^{D} = -\left(t_w^{cem}+t_w^{aem} \right)\frac{i(x)}{F}-\left(L^{cem} \left(p_{s, osm}^{cem, L}(x)-p_{s, osm}^{cem, R}(x) \right)+L^{aem} \left(p_{s, osm}^{aem, R}(x)-p_{s, osm}^{aem, L}(x) \right)\right)\frac{\rho_w}{M_w}$ has_Nernst_diffusion_layer==True

Note

1 ϕm(x), ϕdohm(x) or ϕdnonohm(x) takes 0 if its corresponding configuration is turned off (value == False).

Frictional pressure drop

This model can optionally calculate pressured drops along the flow path in the diluate and concentrate channels through config has_pressure_change and pressure_drop_method. Under the assumption of identical diluate and concentrate channels and starting flow rates, the flow velocities in the two channels are approximated equal and invariant over the channel length when calculating the frictional pressure drops. This approximation is based on the evaluation that the actual velocity variation over the channel length caused by water mass transfer across the consecutive channels leads to negligible errors as compared to the uncertainties carried by the frictional pressure method itself. Table 7 gives essential equations to simulate the pressure drop. Among extensive literatures using these equations, a good reference paper is by Wright et. al., 2018 (References).

Table 7 Essential equations supporting the pressure drop calculation
Description Equation Condition
Frictional pressure drop, Darcy_Weisbach $p_L=f\frac{\rho v^2}{2d_H}$ 1 has_pressure_change == True and pressure_drop_method == PressureDropMethod.Darcy_Weisbach
pL= user-input constant has_pressure_change == True and pressure_drop_method == PressureDropMethod.Experimental
Hydraulic diameter $d_H=\frac{2db\epsilon}{d+b}$ hydraulic_diameter_method == HydraulicDiameterMethod.conventional
$d_H=\frac{4\epsilon}{\frac{2}{h}+(1-\epsilon)S_{v,sp}}$ hydraulic_diameter_method == HydraulicDiameterMethod.spacer_specific_area_known
Renold number $Re=\frac{\rho v d_H}{\mu}$ has_pressure_change == True or limiting_current_density_method == LimitingCurrentDensityMethod.Theoretical
Schmidt number $Sc=\frac{\mu}{\rho D_b}$ has_pressure_change == True or limiting_current_density_method == LimitingCurrentDensityMethod.Theoretical
Sherwood number Sh = 0.29Re0.5Sc0.33 has_pressure_change == True or limiting_current_density_method == LimitingCurrentDensityMethod.Theoretical
Darcy's frictional factor f = 4 × 50.6ϵ − 7.06Re − 1 friction_factor_method == FrictionFactorMethod.Gurreri
f = 4 × 9.6ϵ − 1Re − 0.5 friction_factor_method == FrictionFactorMethod.Kuroda
Pressure balance pin − pLl = pout has_pressure_change == True

Note

1 As discussed in the last paragraph, in this section we assumed a constant linear velocity (in the cell length direction), v, in both channels and along the flow path. This v is calculated from the volume flow rate at x = 0 by the property package.

Nomenclature

Table 8. Nomenclature
Symbol Description Unit
Parameters
ρw Mass density of water kg m − 3
Mw Molecular weight of water kg mol − 1
Variables
N Molar flow rate of a component mol s − 1
J Molar flux of a component mol m − 2s − 1
b Cell/membrane width m
l Cell/membrane length m
t Ion transport number dimensionless
I Current input A
i Current density Am − 2
U Voltage input over a stack V
u x-dependent voltage over a stack V
n Cell pair number dimensionless
ξ Current utilization coefficient (including ion diffusion and water electroosmosis) dimensionless
z Ion charge dimensionless
F Faraday constant C mol − 1
D Ion Diffusivity m2s − 1
δ membrane thickness m
c 1 Solute concentration mol m − 3
tw Water electroosmotic transport number dimensionless
L Water permeability (osmosis) ms − 1Pa − 1
posm Osmotic pressure Pa
rtot Total areal resistance Ωm2
r Membrane areal resistance Ωm2
rel Electrode areal resistance Ωm2
d Channel height m
κ Solution conductivity Sm − 1 or Ω − 1m − 1
η Current efficiency for desalination dimensionless
P Power consumption W
PQ Specific power consumption kW h m − 3
Q Volume flow rate m3s − 1
ϕm Nonohmic potential across a membrane V
ϕdohm Ohmic potential across a Nernst diffusion layer V
ϕdnonohm Nonohmic potential across a Nernst diffusion layer V
Δ Nernst diffusion layer thickness m
Db General mass diffusivity of the salt molecule in the bulk solution m2s − 1
ilim Limiting current density Am − 2
ilim, 0 Limiting current density experimentally evaluated from the untreated salt solution Am − 2
λ Equivalent conductivity of the solution m2Ω − 1mol − 1
A Experimental parameter to calculate ilim Cmol − 1m1 − BsB − 1
B Experimental parameter to calculate ilim dimensionless
dH Hydraulic diameter m
ϵ Spacer porosity dimensionless
f Darcy's fraction factor dimensionless
pL Frictional pressure drop per unit of length Pam − 1
v Flow linear velocity ms − 1
ρ Aqueous phase density kgm − 3
Sv, sp Volume-specific surface area of the spacer m − 1
μ Dynamic viscosity Pas − 1
Re Renold number dimensionless
Sc Schmidt number dimensionless
Sh Sherwood number dimensionless
Subscripts and superscripts
C Concentrate channel
D Diluate channel
j Component index
in Inlet
out Outlet
cem Cation exchange membrane
aem Anion exchange membrane
iem Ion exchange membrane, i.e., cem or aem
L The left side of a membrane, facing the cathode
R The right side of a membrane, facing the anode
s location of the membrane surface
b location of bulk solution
+ or − mono-cation or mono-anion

Note

1 When no component subscript (j), c refers to the concentration of the salt molecules (e.g., NaCl)

References

Strathmann, H. (2010). Electrodialysis, a mature technology with a multitude of new applications. Desalination, 264(3), 268-288.

Strathmann, H. (2004). Ion-exchange membrane separation processes. Elsevier. Ch. 4.

Campione, A., Cipollina, A., Bogle, I. D. L., Gurreri, L., Tamburini, A., Tedesco, M., & Micale, G. (2019). A hierarchical model for novel schemes of electrodialysis desalination. Desalination, 465, 79-93.

Campione, A., Gurreri, L., Ciofalo, M., Micale, G., Tamburini, A., & Cipollina, A. (2018). Electrodialysis for water desalination: A critical assessment of recent developments on process fundamentals, models and applications. Desalination, 434, 121-160.

Spiegler, K. S. (1971). Polarization at ion exchange membrane-solution interfaces. Desalination, 9(4), 367-385.

Wright, N. C., Shah, S. R., & Amrose, S. E. (2018). A robust model of brackish water electrodialysis desalination with experimental comparison at different size scales. Desalination, 443, 27-43.