Analysis of Membrane Transport Equations for Reverse Electrodialysis (RED) Using Irreversible Thermodynamics.

Reverse electrodialysis (RED) is an electro-membrane process for the conversion of mixing energy into electricity. One important problem researchers' face when modeling the RED process is the choice of the proper membrane transport equations. In this study, using experimental data that describe the membrane Nafion 120 in contact with NaCl aqueous solutions, the linear transport equation of irreversible thermodynamics was applied to calculate the power density of the RED system. Various simplifying assumptions about transport equation (i.e., four-, three-, and two-coefficients approaches) are proposed and discussed. We found that the two-coefficients approach, using the membrane conductivity and the apparent transport number of ions, describes the power density with good accuracy. In addition, the influence of the membrane thickness and the concentration polarization on the power density is also demonstrated.

1. Introduction

Global population growth contributes to the escalation of energy demand for households, transportation, and industry. The use of fossil fuels results in various environmental threats. Therefore, the search is increasingly widespread for alternative methods of energy conversion, e.g., osmotic into electrical [1] and mechanical energy (so called “osmotic pump” [2,3], or the conversion of electrokinetic energy [4]. Reverse electrodialysis (RED) is an electro-membrane process for the conversion of mixing energy into electricity. The concept of RED was invented in 1950 [5,6]; however, it was not heavily investigated until the 21st century, when this technique regained the attention of researchers looking for novel, unconventional, and renewable energy sources [7,8,9].

One of the important RED issues is the cost-effectiveness of the process. Turek and Bandura [10] suggested that energy conversion using RED is not profitable; however, Giacalone et al. [11] presented a much more optimistic perspective. There are a number of parameters that influence the efficiency of RED, like ionic shortcut currents [12], type of electrodes [13], kind of membrane [14], membrane thickness [15], fouling issues [16,17], efficiency of solutions mixing in the stack [18,19], and geometry of spacers [20].

In this section, we present a description of membrane transport in RED. In the literature discussing the modeling of RED, apart from very simplified membrane transport (i.e., considering only the transport numbers of ions [21]), usually the Nernst–Planck equation or its transformations are used [22,23,24,25,26,27,28]. In some research, the volume flow caused by the osmotic and electroosmotic transport of water through the membrane was considered [25,29]. Whereas in other research, only osmotic transport was considered [22,24]. In other studies, no volume transport was anticipated [27]. The problem of membrane transport is probably neglected because of the many other important aspects of RED. A general (i.e., not exactly related to RED) review of transport through ion-exchange membranes can be found elsewhere [30].

The Nernst–Planck equation is a simplified version of the linear transport equation of irreversible thermodynamics. The transport equation is comprised of six transport coefficients for the system membrane in contact with a single electrolyte aqueous solution, under the condition that the Onsager symmetry of the coefficients is fulfilled [31]. The determination of all these coefficients requires extensive work especially because their concentration dependence should be considered [32,33,34,35]. Consequently, the practical use of the full transport equation is limited.

Therefore, the aim of this work was to discuss which of the simplifications of the irreversible thermodynamics transport equation (ITE) is useful for calculating the power density in RED. The following three approaches based on four, three, and two coefficients (to be determined experimentally) were analyzed and discussed in detail:

the four-coefficients approach—the pressure difference in the full set of ITE is neglected,

the three-coefficients approach—the coupling of ion flows with the gradient of chemical potential of water is neglected,

the two-coefficients approach (i.e., Nernst–Planck equation)—the ion flux is a function of the gradient of its electrochemical potential only.

The results of these simplifications (i.e., values of power density vs. electric current density) are compared with the values obtained from the full set of ITE. The transport coefficients used in these latter equations were determined for the cation-exchange membrane Nafion 120 in contact with NaCl solutions [34,35] in a wide range of concentrations (0.05–3.7 M). To consider the concentration dependence, the differential form of the transport equation was used. This method was applied for the first time by McCallum and Meares to describe transport through an ion-exchange membrane [36]. The Nafion 120 membrane is relatively thick (0.3 mm); therefore, the transport coefficients are rescaled to smaller thickness values to demonstrate the influence of that parameter on the power density.

2. Theory

In the following, one basic part of the RED module is considered, i.e., one ion-exchange membrane with adjacent concentration polarization layers (Figure 1). The steady-state conditions in the system are assumed, i.e., the fluxes of ions in the membrane and both polarization layers are the same.

Figure 1

The system: ion-exchange membrane with adjacent boundary diffusion layers. lp—thickness of concentration polarization layer; lm—thickness of the membrane; c′, c″—concentrations of ions in the bulk solution; c′(m), c″(m)—concentrations of ions at the membrane surface.

2.1. Transport Equations of Irreversible Thermodynamics (ITE)

The starting point for the ITE approach to isothermal conditions is the dissipation function, Φ, the product of entropy production, and temperature, which can be expressed as the sum of driving forces, X, and coupled with fluxes, J:

For discontinuous systems, the forces are represented by the differences in chemical potentials of transported species. In continuous systems, the negative values of gradients of these potentials are used.

According to irreversible thermodynamics, it is assumed that fluxes and forces are linearly related: and the Onsager reciprocity relation L = L holds if the forces are not too high.

One can choose many sets of fluxes and forces; all of them have to fulfill the condition of the entropy production invariance:

From this condition, we obtain various forms of the fluxes and forces, among which is the set of practical equations useful for experimental determination of transport coefficients [37]: where J1 is the flux of ions 1 (cations), ν1 is the number of ions 1 in the electrolyte molecule, s is volume flux, I is electric current density, Δp is the pressure difference, Δπ is the osmotic pressure difference, is the mean concentration of electrolyte solutions contacting the membrane (c′, c″), and is the mean transport coefficient for that concentration range. ΔE is the electric potential difference measured with electrodes reversible to anions (2):

The power density is calculated from the electric potential difference Δφ, determined from Equation (4) assuming that .

The determination of all L coefficients is a laborious task. In Table 1, we list the experiments we performed to determine L in Equation (3).

Table 1

List of experiments, equation for 1:1 electrolyte.

Experiment	Equation
conductivity	κ≡l_mIΔE_c,p=l_mL_EE
the real transport number of ion 1	t¯_1≡FJ_1I_Δc,Δp=0=L_πEL_EE
electroosmotic coefficient	W≡FJ_vI_Δc,Δp=0=v¯_0t¯_0+v¯_st¯_1=L_pEL_EE
the apparent transport number of ion 1 (electromotive force)	t¯_1,app=FΔEΔμ_s_I,Δp=0=F1−c˜_sv¯_sL_EπL_EE−c˜_sL_EpL_EE t¯_1,app_Δc→0=t¯_1−c_sc_0t¯_0
electrolyte permeability	L_s≡J_sΔπ_I,Δp=0=L_ππ/c˜_s−L_πp+L_EpL_πE−L_πEL_Eπ/c˜_sL_EE
osmotic permeability	L_v,osm≡−J_vΔπ_I,Δp=0=L_pp−L_πp/c˜_s+L_pEL_Eπ/c˜_s−L_pEL_EpL_EE
hydrodynamic permeability	L_p≡J_vΔp_I,Δc=0=L_pp−L_pEL_EpL_EE

Having these coefficients determined for the vanishingly small concentration differences (Δπ, Δμ → 0) allowed us to assume the Onsager reciprocity relation (L = L) and determine the concentration dependence of six L (α, β = π, p, E) [34,35]. The method of such determination was described in detail elsewhere [32]. The finite concentration differences (L) in Equation (3) are not symmetric because of their strong concentration dependence [32,38].

Calculating the ion flux, volume flux, and electric potential difference, considering the concentration dependence of L is the most exact approach. To calculate these quantities for the assumed conditions (concentrations on both sides of membrane, electric current density, membrane thickness), we applied the method proposed by McCallum and Meares [36]. According to this method, the membrane is divided into nm slices and for each slice the transport is given by the set of Equation (3) with local variables (electrolyte 1:1): where: where l is the thickness of the real membrane for which the coefficients L = f(c) were determined and w is the width of slice equal to the assumed membrane thickness (l,assumed) divided by the number of slices (w = l,assumed/nm). The concentration difference is expressed as Δπ and the formula .

The fluxes across each slice are the same, assuming the stationary state in the considered system. Therefore, for nm slices, we have 3n equations and 3nm unknowns: c, p, E, k = 1, …, n − 1, J1, and J, which we solve for given values of parameters (electric current, concentration, pressure on both sides of membrane, and the assumed membrane thickness). All of these refer to the hypothetical (virtual) solution being in equilibrium with the membrane at coordinate x = kw, k = 1, …, n − 1; such a solution is defined by equality of chemical potentials of the transported species in that solution and in the membrane at a given x.

The concentration polarization layers adjacent to membrane surfaces were treated in a similar way as the membrane. The experimental data for NaCl solutions were taken from [39]. The used equations are provided in Appendix A.

A comment should be made on the extended Nernst–Planck equation, frequently used in the literature for describing the membrane transport of ions. Written in terms of hypothetical solution parameters, it takes the following form: where K is the partition coefficient of ion I and α is the convective coupling coefficient [40,41]. There are 2 coefficients per ion (K is treated as one coefficient); thus, for two ions and J described by 3 coefficients (L, L, L), there are 7 coefficients—more than in the set of Equation (3) assuming the symmetry of coefficients. For this reason, the extended Nernst–Planck equation will not be discussed here.

2.2. The Four-Coefficients Approach

In this case, the pressure difference across each slice is assumed to be zero. Although on both sides of membrane the pressure is the same, it may change across the membrane to maintain the volume flux across each slice the same. As the number of unknowns is reduced, the equation for J in Equation (5) can be neglected (J1 is needed in the concentration polarization layers). Thus Equation (5) is reduced to: where , .

The four coefficients (L′, L, L′, L) can be determined from P, , k, and —the Electromotive Force (EMF) measurements. However, practically, it is better to determine the electroosmotic coefficient W (Table 1) and to calculate from W and .

2.3. The Three-Coefficients Approach

In the equation of ion transport, the coupling with water, i.e., the term L0X0 is neglected. Thus, Equation (1) with L replaced by P/RT reduces to: where X is:

We assume that L21 = L12; therefore, there are 3 parameters that can be obtained from (1) specific conductivity, km, (2) electrolyte permeability P, and (3) the apparent transport number, . Instead of , we could choose the real transport number of cation, ; however, because includes the electroosmotic number of water and is much easier to determine than , we choose . Notably, and obtained from Equation (7) are the same because the flux of water is neglected here (see formulae for and in Table 1). The EMF measurements allow for the determination of the apparent transport number of ion 1. Using Equation (7), we obtain the following set of practical equations (for the 1:1 electrolyte): where P1 is given by: and P by:

2.4. The Two-Coefficients Approach

In this approach, we assume that the flux of ion I depends only on the gradient of electrochemical potential of that ion:

The following practical equations are derived from Equation (10):

There are 2 transport coefficients only, P1 and P2, with four experimentally available parameters: (1) specific conductivity k, (2) transport number of cation , (3) electrolyte permeability P, and (4) . The relation between these quantities and P1, P2 is as follows (assuming 1:1 electrolyte, and denoting ion 1 as a cation): and

Similar to the 3-coefficients approach, the expressions for and are the same because the water transport is also neglected.

The determination of by the EMF method is much easier and therefore much more popular than the determination of the real transport number, . We chose the following pairs of experiments to determine P1 and P2:

(1) k, , (2) k, P, (3) , P.

Additionally, contrary to , the value of includes the electroosmotic transport of water; thus, to some extent, the transport of water molecules is considered.

3. Results and Discussion

As mentioned above, in this work, only one membrane (cation-exchange membrane) with adjacent concentration polarization layers is discussed.

In Figure 2, a typical dependence of the power density P on the current density I, calculated as P = −IΔφ, based on the full transport equation (Equation (5)), is shown for two concentration differences (0.1–1.0 M and 0.05–3.7 M) and two values of membrane thickness (0.1 and 0.3 mm). Two common trends were observed: the higher the concentrations difference, the higher the value of P; the lower the membrane thickness, the higher the value of P.

Figure 2

Electric power vs. −I calculated using Equation (5) (full ITE) for the real thickness of Nafion 120 equal to l = 0.3 mm, and the assumed value l = 0.1 mm, for 0.1–1.0 M and 0.05–3.7 M NaCl, T = 298 K. l—thickness of the concentration polarization layer, l—thickness of the membrane.

The four-coefficients approach, neglecting the pressure changes (Equation (6)), yields similar results. For 0.1–1.0 M NaCl, the difference in Pmax is below 1% irrespective of the membrane thickness. For higher concentration differences (0.05–3.7 M), this difference increases to ca. 5–6%.

Changes of pressure in the hypothetical/virtual solution affect the ion fluxes. Some examples of pressure variations are shown in Figure 3. The negative values of pressure indicated that there must have been an underpressure inside the membrane needed to keep the volume flow constant. Similar results were obtained by McCallum and Meares [36]. They suggested that the membrane thickness could decrease in that case. At higher I values, the pressure becomes positive. This can be explained by different changes in the profiles of electroosmotic and osmotic parts of J across the membrane caused by the change of concentration profiles (Figure 4). To keep J constant, the pressure part of J should change accordingly, which demands various profiles of pressure. Assuming the constant pressure value (the four-coefficients approach), the volume flux across the membrane is not constant—it can deviate from the actual steady state value even several times, which is unrealistic.

Figure 3

Pressure profile across the membrane calculated from Equation (5) (full ITE) for l = 0.3 mm and l = 0.1 mm, for 0.1–1.0 M and 0.05–3.7 M NaCl, T = 298 K. Thickness of the concentration polarization layer l = 0.05 mm. Thickness of the membrane—l.

Figure 4

Profiles of electroosmotic (el), osmotic (osm), and pressure-driven (p), components of J calculated from Equation (5) (full ITE) for the current density –I = 100 and 200 A/m2; one can notice a dominant role of electroosmotic contribution to J at higher current density.

The influence of the membrane thickness on the maximum of P and the current density at Pmax is shown in Figure 5. On both curves (Pmax = f(l) and −I(Pmax) = f(l)), maxima are observed at l < 0.1 mm, and the maximum on Pmax is at a slightly higher l value than that of −I(Pmax). The higher concentration difference demands a thicker membrane (maxima are shifted toward thicker membranes). The thinner the membrane, the higher I, but also the higher the permeability of electrolytes; the impact of diffusion layers (of the assumed constant thickness) also becomes relatively stronger.

Figure 5

Dependence of the maximum of power density, Pmax, and the current density, −I, at Pmax on the membrane thickness, l, for 0.1–1.0 M and 0.05–3.7 M NaCl, T = 298 K, Equation (5), full ITE.

The dependence of the maximum of power density Pmax and the current density at Pmax on the diffusion boundary layer thickness l is a monotonically decreasing function (Figure 6). Due to a substantial decrease in Pmax with l, the reduction of the thickness of that layer by ensuring adequate stirring of the solutions at the membrane surface is an important practical issue. Even more important is the decrease in l on the lower concentration side of membrane. To reduce the concentration polarization effect, an extra energy for pumping is needed.

Figure 6

Dependence of the maximum of power density, Pmax, and the current density, −I, at Pmax on the diffusion boundary layer thickness, l, for 0.1–1.0 M and 0.05–3.7 M NaCl, T = 298 K, Equation (5), full ITE.

The Two-Coefficients Approach

The results for the two-coefficients approach, for which P1 and P2 are calculated either from and k (TC), or from and P (TP), or from km and P (CP), are shown in Figure 7. Among these three combinations of coefficients, the smallest deviation from the power density calculated from the full ITE was observed for the TC combination (, k). It yielded even better results than the three-coefficients approach based on all three parameters (TCP: , k and P).

Figure 7

Power density vs. current density for full model and approximate approaches: TCP—3-coefficients approach (P1, P2, and L12 calculated from , k, and Ps), (t1, P); TP, CP, TC—2-coefficients approach (P1 and P2 calculated from and Ps, k and P, and k, respectively); the Nafion 120 membrane separates (a) 0.1 and 1 M NaCl, (b) 0.05 and 3.7 M NaCl; l = 0.1 mm, l = 0.05 mm, T = 298 K.

Because of the quite good fit of the TC approach, another question arose: is it necessary to determine the concentration dependence of and k in the whole range c′–c″, or is it enough to determine in that range and km for the mean concentration of c′ and c″?

In that case (P1, P2 = const), we can determine the integrated form for J1 from Equation (10):

ΔE can be expressed as a function of a certain mean concentration of c′ and c″:

P1 and P2 are determined from km (Equation (12)), with cs replaced by , and calculated from (15) with dE/dμs replaced by (E(c″) − E(c′))/(μs(c″) − μs(c′)). To check this simplification, the calculations were performed for being the arithmetic and logarithmic means of c′, c″. The boundary diffusion layers were treated as in the previous cases to keep the same conditions for all discussed approaches. To simplify calculations, in the calculations of P1, P2, and , the external bulk concentrations were used, not those at the membrane surfaces. A comparison of this approach (TC const.) with the full ITE equation and with the TC approach is shown in Figure 8.

Figure 8

Power density vs. current density for full model—comparison of the TC const approach with logarithmic mean with the TC approach and the full ITE for the concentrations 0.1–1 M and 0.05–3.7 M NaCl; l = 0.1 mm, l = 0.05 mm, T = 298 K.

This simplification of the logarithmic mean of concentration (marked in Figure 8 as “TC const ln”) yields comparable results to the other approaches for the lower concentration difference (0.1–1.0 M). For 0.05–3.7 M the predictions of “TC const ln” are too low at the current density range 0–j(Pmax).

In this case, the concentration dependence of and k is needed to improve the predictions of power density. The arithmetic mean of concentration yields much worse results especially for high concentration difference and is not shown in Figure 8.

4. Conclusions

Among the tested simplifications of the full set of transport equations of irreversible thermodynamics (six independent coefficients), the four-coefficients approach (k, , , P,dif) with only the pressure effect neglected yields almost the same results for 0.1–1 M as the full ITE approach.

The three-coefficients approach (k, , P,dif), where only ion fluxes are considered with their couplings, slightly overestimates the power density at lower concentration difference. However, at higher concentration differences, it strongly overestimates and becomes much less accurate than the two-coefficients approach (equation of the Nernst–Planck type, assuming no coupling between ion fluxes) based on k and . Other combinations of the two-coefficients approach are substantially worse. The model based on k and can be further simplified by using constants P1 and P2 calculated from the mean values of k and in a given concentration range c′–c″ and using the logarithmic mean of these concentrations. These coefficients are easily accessible experimentally.

None of the discussed approximations considered the volume flow effect on the concentrations of the solutions contacting the membrane. To do so, we need to estimate the osmotic volume permeability and the electroosmotic coefficient W.

To obtain the optimal power density for given concentrations of feed solutions, the appropriate membrane thickness should be chosen. The higher the concentration difference, the thicker the membrane should be. As the position of Pmax depends on the membrane properties, the optimal thickness should be independently determined for cation- and anion-exchange membranes.

The thickness of the diffusion boundary layer should be as low as possible especially on the low concentration side of the membrane—it strongly decreases the power density.