Aspects of Mathematical Modelling of Pressure Retarded Osmosis.

In power generating terms, a pressure retarded osmosis (PRO) energy generating plant, on a river entering a sea or ocean, is equivalent to a hydroelectric dam with a height of about 60 meters. Therefore, PRO can add significantly to existing renewable power generation capacity if economical constrains of the method are resolved. PRO energy generation relies on a semipermeable membrane that is permeable to water and impermeable to salt. Mathematical modelling plays an important part in understanding flows of water and salt near and across semipermeable membranes and helps to optimize PRO energy generation. Therefore, the modelling can help realizing PRO energy generation potential. In this work, a few aspects of mathematical modelling of the PRO process are reviewed and discussed.

1. Introduction

When fresh river water mixes with salty sea water, energy equal to the one that is produced in a waterfall about 200 meters high is lost [1]. If only 25%–30% of this energy is converted to electricity it will be equivalent to damming a river with a hydroelectric dam with a height of about 60 meters. The concept of using salinity gradient for power generation was first introduced by Pattle [1] in 1954 who described “hydroelectric pile” apparatus (now called reversed electrodialysis) to harness this energy. Hitherto, it appears that a method based on osmosis, proposed by Loeb in 1973 and first published in 1975 [2], is closer to practical realization. The first, and hitherto the only power plant prototype based on pressure retarded osmosis (PRO) was commissioned by Norwegian state-owned power company Statkraft in 2009 and is now in operation.

Generating energy based on PRO relies on creating a pressure in a more salty draw solution using the osmotic flow of water from a less salty feed solution through a semipermeable membrane (see Figure 1). Energy is produced by a generator using the “excess” pressure and flow due to the osmotic flow of water.

Figure 1

Schematic representation of a classic pressure retarded osmosis (PRO) energy generation plant with a compartmental geometry (adapted from [3]). The feed and draw compartments are assumed to be well stirred, so that the concentration at the membrane surface is the same as the outflowing concentration: and .

PRO energy generation has recently been extended to include closed loop systems, for example, one based on ammonia-carbon dioxide [4] and novel PRO coupled with electric power generation from the electrokinetic streaming potential [5]. A PRO system capable of utilizing low-grade heat not recoverable with existing technologies was presented and investigated [6]. For up to date information on these new developments as well as the applicability and practical viability and current limitations of PRO systems readers are referred to recent comprehensive reviews [7,8,9,10,11].

This short review aims to describe mathematical modelling most relevant to energy generation using PRO. Mathematical modelling where analytical approach is possible is given a priority with mathematical derivations somewhat more detailed than usual in the PRO literature. It is hoped that this detailed description of mathematical solutions will be beneficial for students and scientists new to this area of research who need to understand the details of mathematical modelling.

2. Ideal Membrane

2.1. Compartmental Configuration

The power output of the PRO power generation scheme presented in Figure 1 is the difference between the power produced by the generator (W) and that used by the pump (W). Assuming that the pressure difference across the pump and the generator (, also called working pressure) is the same and is equal to the hydrostatic pressure difference across the membrane, the generator and pump powers can be expressed as: and respectively, where is the osmotic water flux, is the surface area of the membrane and is the draw solution pump flow rate or flow rate into the draw compartment. Therefore, the power generated per unit area of the membrane () is [3]:

The osmotic water flux is positive (that is, from the feed to the draw solution) if the osmotic pressure difference across the membrane () is greater than the hydrostatic pressure () and is linearly proportional to the difference between the osmotic and hydrostatic pressures [3]: where is the water permeation coefficient of the semipermeable membrane. It is important to note that Equation (1) does not account for losses due to the low pressure water pumping (feed solution side) required for the operation of a PRO plant and inefficiencies of converting energy at the generator and the high pressure pump (draw solution side).

Combining Equations (1) and (2) yields:

Since the hydrostatic pressure () can be varied, W can be maximized with regard to . Differentiating Equation (3) with respect to and finding where the derivative of power is zero () gives optimal hydrostatic pressure: . Therefore, the maximum power is:

Equation (4) indicates ways to increase power output of PRO, which is possible by selecting membranes with higher permeability (), or increasing osmotic pressure difference (), for example, by increasing salt concentration of the draw solution. Increasing the osmotic pressure difference is especially attractive as dependence of the maximum power in (4) on is quadratic, but this option can be limited by what solutions are available.

The osmotic pressure difference is determined by the concentration difference across the membrane: where and are concentrations of salt on the feed and draw sides of the membrane respectively and is a monotonically increasing function (see Equation (7)). For the compartmental configuration it is assumed that, due to the effective stirring, concentrations near membrane surfaces are the same as bulk concentrations, which in turn are equal to the outflowing concentrations: and . Due to the osmotic water flow through the membrane (), these concentrations are not equal to the feed and draw concentrations flawing into the compartments () and () respectively, but can be related to them by considering mass conservation of water ( and ) and salt ( and ) in the system, yielding: where and are the feed solution flow rates into and out of the feed compartment and it is assumed that the draw and feed compartments are well stirred and there is no concentration polarization (see Section 3) at the membrane surface. Increasing the difference between concentrations of the feed and draw sides of the membrane () increases and therefore . The power can be increased further by increasing the feed and draw pump flow rates, as according to Equation (6) this will increase . The latter could be counterproductive if losses caused by the increased pumping rates are higher than gains due to the rise in the osmotic pressure.

The osmotic pressure is linearly proportional to a salt concentration (C, expressed in moles per volume) when the concentration is not large [12]: where R is the ideal gas constant, T is the absolute temperature and i is the dimensionless van't Hoff factor. The van’t Hoff factor for electrolytes depends on the degree of the dissociation and the number of ions. For NaCl solutions of up to the concertation of sea water i ≈1.9 was suggested [12].

If the feed and draw concentrations are known, Equations (4) to (7) can be used to calculate the maximum power per unit area of a membrane for a given membrane permeability.

2.2. Counterflow Configuration

Compartmental geometry is easy to realize and analyze in laboratory conditions and where only a small area of the membrane is necessary. In this case, draw and feed compartments volumes to the surface of the membrane ratio are not small and it is easy to arrange the effective stirring so that the compartments are well stirred. When stirring is effective, the concertation inside the compartments is uniform and practically the same near the input and the output of the draw and feed solutions (top and bottom in Figure 1 respectively). In practical applications of PRO very large areas of membrane are required: millions of square meters of the membrane to generate just tens of megawatts of power. In this case, the effective stirring is not practical and a concentration gradient will form between the input and the output of the draw and feed solutions due to the water flow through the membrane. Schematic diagram of a PRO unit with the counterflow configuration is presented in Figure 2. The tubular membrane shape is perhaps the most practical similar to desalination units where membranes are organized as multiple small bore tubes that are placed in a much larger metal tube casing. The small diameter of the tubes leads to a large surface area of the membrane in a relatively small volume.

Figure 2

(a) A diagram representing tubular geometry and counterflow configuration for PRO. (b) Simplified representation of the tubular geometry as a square channel. Concentrations ( and ) and flows (, and ) are the functions of the position x along the channel for the counterflow configuration.

In Figure 2b the tubular geometry is replaced with the square channel representation for simplicity and better visualization of the problem. To derive transport equations for this counterflow configuration we will generally follow the approach presented by Sharqawy et al. [13], changing to notations consistent with this work. As concentrations and flows in the feed and draw channels are now functions of x, where x is the distance along the channel, Equation (2) for the water flow across the membrane becomes: where , and are the draw and feed concentrations and the water flow per unit area of the membrane at position x along the channel respectively. Here we assumed that the hydrostatic pressure, , is constant along the channel, even though pressure gradient may form in a thin tube due to the water viscosity, this pressure is likely to be small compared to the hydrostatic pressure. We have also assumed in Equation (8) that concentrations in the feed and draw channels are uniform across the channel. This is only possible if there is no concentration polarization (see Section 3).

We further assume, as for the compartmental geometry, that only the pure water flows through the membrane. Taking into account the pure water flow through the section of the membrane from x to away from the feed to the draw solution the feed flow rate can be expressed as: where is the area of the membrane between x and x + dx. This area can be further expressed as , where w is the width of the channel or the circumference of the tubular membrane.

Substituting in Equation (9) yields:

Using the initial condition and integrating Equation (10) gives: where the integral in Equation (11) was replaced with , the flow of pure water through the membrane from feed to draw solution between 0 and x.

A similar consideration for the draw flow rate leads to the differential equation identical to Equation (10), with a change of . Using the initial condition , where L is the length of the channel, yields:

Here we assumed for simplicity that is positive, even though is directed against the positive direction of the x axis (Figure 2b).

As there is no salt flow through the membrane the flux of salt for the feed and draw sides is constant: and . Using these equations and Equations (11) and (12), concentrations in the feed and draw solutions can be expressed as:

In Equations (13) parameters and are operational parameters of the PRO system that are determined, and the only undetermined function is . The differential equation for this function can be derived by noting that and using Equations (8) and (13):

Equation (14) is a separable first order differential equation which can be solved by integration:

The integral in this equation can only be integrated analytically when the osmotic pressure is linearly proportional to a salt concentration (see Equation (7)). In this case the integrand in Equation (15) can be presented as two simple fractions of the form , where a and b are constants, and the integration yields logarithmic expressions (for details see Equation [13]). Even in this simpler linear case it is not possible to represent explicitly and numerical approach is required. Sharqawy et al. [13] analyzed in detail the linear osmotic pressure case by dimensionalising Equation (15) and introducing notations analogous to a mathematical modelling of heat exchangers. The non-linear osmotic pressure case was also considered for both the counterflow and parallel flow ( and in the same direction) configurations, but it was found that for the case of the seawater as the draw solution and the river water as the feed solution the error of using the linear approximation is not significant [13]. Sharqawy et al. [13] also modelled parallel-flow configuration for PRO, but concluded that, just as for the heat exchanges, that the counterflow configuration is more efficient for the same area of the membrane used. Approach similar to work of Sharqawy et al. [13] was used to analyze limits of power generation due to finite membrane area [14].

After numerically solving Equation (15), the total power generated by the PRO system () can be determined as: and further analyzed to maximize the power relative to the operational parameters of the PRO system, including the hydrostatic pressure (). Given the number of parameters, this is not a straightforward exercise and requires a numerical approach.

3. Concentration Polarization

Concentration polarization is arguably the most significant problem that dramatically reduces the power output of the PRO process and reduces the applicability of equations presented above for the ideal membrane. The polarization had been experimentally investigated and mathematically modelled for the first time more than 30 years ago [3,15]. Experimental work on the polarization was later conducted with more modern membranes [16]. In this review the mathematical modelling of the polarization will be presented for a realistic asymmetric membrane with the porous support layer against the feed solution as shown in Figure 3. The diagram of the concentration polarization is shown in Figure 3 and is represented by the significant increase in the feed solution concentration at the surface of the active layer. This increase is due to the convective transport of salt by the water flow in the support layer of membrane, as the salt cannot penetrate through the active layer of membrane and needs to diffuse against the flow. It is also shown that concentration changes in the unstirred boundary layers near the membrane surfaces (Figure 3). Mathematical modelling of concentration polarization that takes into account the unstirred boundary layers was presented in [17]. The concentration polarization in this case is inversely proportional to the Sherwood number, which in turn depends on the Reynolds number, Schmidt number and geometric dimensions of the channel [18]. Salt leakage/flow across the active layer of the membrane in the direction opposite to the water flow due to the membrane imperfections further increases the concentration polarization.

Figure 3

A diagram representing the concentration polarization for an asymmetric membrane with the support layer due to the water flow through the membrane. Here C and C are concentrations of salt on the feed and draw sides of the active layer respectively and C and C are the bulk concentrations on the feed and draw sides of the membrane respectively.

For simplicity, let’s assume that the concentration changes in the unstirred boundary layers near the membrane surfaces are negligible. Then the boundary condition at the support layer will be: [3] where it was assumed that x = 0 at the left boundary of the support layer and the positive x is in the direction of the osmotic water flow.

The salt flux in the support layer is the sum of diffusive and convective fluxes: [17] where D is the diffusion coefficient of salt and ε is the porosity of the support layer. For the steady state is constant, and Equation (18) is the first order inhomogeneous differential equation with constant coefficients and can be easily integrated yielding: where E is the integrating constant which can be determined using boundary condition (Equation (17)), so that:

The concentration at the feed side of the active layer can now be determined as , where t is the thickness and τ the tortuosity of the support layer, therefore: where is a measure of diffusional resistance to salt transport in the support layer. The support layer structural parameter S () depends only on the microstructure of the support membrane and is now commonly used in the PRO literature [17,19,20] to describe the concentration polarization. It follows from definition of these parameter that K = S/D.

It is appropriate to assume that the salt leakage across the active layer () is linearly proportional to the salt concentration difference across the layer: [3] where B is the salt permeation coefficient across the active layer. Substituting from Equations (22) to (21) yields linear equation for that can be easily solved:

In this equation, the concentration at the feed side of the active layer is expressed using parameters of the membrane, the bulk salt concentrations and osmotic water flow (). Equation (23) together with Equations (2) and (5) form an explicit set of equations which allow determination of from parameters of the membrane, the bulk salt concentrations and the hydrostatic pressure (). These equations generally need to be solved numerically.

Lee et al. [3] using equations similar to Equation (23), Equation (2) and Equation (5), assuming the linear relationship between concentration and osmotic pressure (Equation (7)) with pure water as a feed solution () and zero hydrostatic pressure () derived expression for the diffusional resistance (K): where is experimentally measured water flow for and .

The mathematical approach is very similar when the concentration changes in the unstirred boundary layers near the membrane surfaces are not negligible. Readers are referred to works by Yip et al. [17] and McCutcheon et al. [18,21] for equations in this case.

4. Module-Scale Analysis of PRO

In Section 2.1 the optimal hydrostatic pressure for an ideal membrane with simple compartmental geometry was determined. For realistic osmotic membranes with significant concentration polarization the process of finding the optimal hydrostatic pressure in general involves a numerical solution with iterations [22,23]. The problem becomes even more challenging for module-scale analysis of PRO when concentration polarization has to be combined with the realistic flow and membrane configurations. Numerical modelling of parallel-flow configuration for PRO was given for the case of concentration polarization [24]. In this work authors introduced a two-dimensional model similar to that presented for crossflow microfiltration and ultrafiltration [25].

Module-scale analysis of PRO for counterflow (or counter-current) and parallel-flow (co-current) configurations were recently investigated taking into account the unfavorable effects of reverse salt flux, internal concentration polarization, and external concentration polarization [20]. In this case, the water and salt fluxes can be represented as: [17,20] where k (=D/δ, δ – effective unstirred layer thickness) is the mass transfer coefficient of the draw solution and are osmotic pressures corresponding to draw and feed solutions respectively. Yip et al. [17] presented detailed derivations of Equations (25) and (26). They also demonstrated that membrane that balances permeability and selectivity allows to achieve the highest potential peak power density for given feed and draw solution concentrations [17].

These equations have to be solved together with the mass transfer equations, which in the case of the counter-current flow operation are: [20] where s is the relative position along the module represented as membrane area from the draw solution entrance to the position in the module and normalized to the total area of the membrane (S). Note that Equation (28) is similar to Equation (10), but has the positive sign in front of J, as the direction of x selected for Equation (10) is opposite to that for s in Equation (28). The boundary conditions for Equations (27)–(30) are , , and [20].

Straub et al. [20] numerically solved Equations (25)–(30) and analyzed power density (, the power generated per membrane area) and the specific energy (, the energy extracted per total volume of the feed and draw solutions combined) using simplifying assumptions of no pressure loss due to pumping of solutions through the feed and draw channels and no inefficiencies in the pressure exchanger or turbine. The approach allows optimizing operating conditions of a realistic PRO system. It also allowed to determine that the maximum specific energy for the current commercial membranes is 1.1 W/m2, only 15% of the power density available for the small scale compartmental (coupon scale) PRO system [20] and well short of 5 W/m2 necessary to produce osmotic power on commercial basis [26].

5. Other Aspects of Mathematical Modelling of PRO

Another potentially practical configuration for PRO energy system is a spiral wound module [27,28]. The experimental and mathematical modelling research for this configuration was conducted by Xu et al. [28]. The dilution of the bulk draw solution for this case is to some extent similar to the counterflow configuration but is further complicated by two different flow paths, axial and spiral [27].

An important aspect of PRO analysis not covered in this review is the thermodynamic efficiency of a PRO process in terms of energy extraction of the Gibbs free energy of mixing. Readers are referred to works by the Elimelech group [29,30] for comprehensive analyses.

Another important factor in realistic module-scale PRO systems are feed and draw channels’ geometry and a pressure loss due to viscose flow of solutions through the channels. Seppälä and Lampinen [31] derived transport equation for osmosis inside a hollow cylindrical fiber, taking into account the cylindrical geometry and hydrostatic pressure drop in the fiber. They solved the equations numerically and found the optimal values of the initial hydrostatic pressure difference between the feed and draw sides of the fiber [31].

In this review, only concentration of a single solute was considered. Such modelling only fully applies for non-dissociating solutes like glucose. In reality PRO systems are most likely to be based on sea and river waters [2] which have multi-ionic composition. In most cases considering NaCl is sufficient to model operation of the PRO systems [24,32], but as NaCl dissociates in water to Na+ and Cl− ions, that creates two ionic species which can differ in their permeability through the membrane. In general, solution diffusion and electro-migration have to be taken into account for multi-ionic systems [33]. Yaroshchuk et al. [33] considered mathematical modelling of transport of multiple ions where diffusion was coupled to electro-migration and concluded that such modelling was important to understand phenomenon such as negative rejection for some ions in particular that spontaneously arising electric fields may yield much higher NaCl rejection, which could be relevant for energy generating PRO systems.

It is clear from discussions in this work that the mathematical problem of optimizing realistic PRO systems is quite challenging. To help with this optimization Sivertsen et al. [12] introduced Iso-watt diagrams that are relatively easy to understand and use, and allow an evaluation of power per unit area of the membrane on basis of membrane characteristics. These diagrams are generated for realistic membranes with concentration polarization and are a useful tool in optimizing PRO energy generating systems, but so far have not been developed to apply to cases of counterflow configurations.

6. Conclusions

Mathematical modelling plays an important part in experimental analysis, development and optimization of PRO energy generating systems. In this work the mathematical modelling for compartmental and counterflow configurations were reviewed and presented in some detail for the simple case of an ideal membrane without concentration polarization. The equations presented in this work for these configurations can be used for getting some insight into optimizing energy generation through varying parameters of the system, but are limited to an ideal membrane. Operating realistic PRO systems leads to a very significant concentration polarization, especially in the support layer of a membrane. Basic approach to mathematical modelling of the concentration polarization and main concepts were reviewed. Although there are approaches that allow optimization of parameters of a PRO system with the concentration polarization for the simple compartmental configuration, more modelling work is required to consider the more practical counterflow configuration and concentration polarization.