Solar-Driven Humidification–Dehumidification Process for Water

Jul 10, 2019 - Read OnlinePDF (4 MB). Supporting Info (1)»Supporting Information ... Abstract Image. This study investigates a desalination system ba...
0 downloads 0 Views 9MB Size
Subscriber access provided by UNIV OF SOUTHERN INDIANA

Process Systems Engineering

A solar-driven humidification-dehumidification process for water desalination analyzed and optimized via equilibrium theory Paolo Gabrielli, and Marco Mazzotti Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.9b02823 • Publication Date (Web): 10 Jul 2019 Downloaded from pubs.acs.org on July 18, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

A solar-driven humidification-dehumidification process for water desalination analyzed and optimized via equilibrium theory Paolo Gabrielli, Marco Mazzotti



Institute of Process Engineering, ETH Zurich, 8092 Zurich, Switzerland

Abstract This study investigates a desalination system based on a humidification-dehumidification process driven by solar modules and proposes a solution within the framework of equilibrium theory. While the optimal operation of the considered humidification-dehumidification system has been investigated in the past, and despite the simplifications and assumptions characterizing equilibrium theory, the proposed solution provides a thorough understanding of the system operation, which has not been obtained yet. First, it allows explaining the asymmetry between the humidifier and the dehumidifier, which does not result from mass transport limitations, but from the functional dependence of the enthalpy of moist air on temperature, which in turn depends on the functional dependence of the saturation pressure of air. On the one hand, humidifying air becomes easier when increasing the temperature, leading to simple wave interactions, to the corresponding smooth transitions, and to a greater entropy generation in the humidifier. On the other hand, dehumidifying air becomes more difficult when decreasing the temperature, leading to shock interactions, to the corresponding sharp transitions, and to a smaller entropy generation in the dehumidifier. Then, it allows describing the behavior of the overall system, i.e. considering the solar modules connecting dehumidifier and humidifier. The optimal system operation is defined by determining the air-to-water mass flow rate ratio that maximizes the productivity of fresh water for any value of ambient conditions and saline water inlet temperature. Such a maximum productivity curve separate a low-productivity and a highproductivity region, with a discontinuity causing a sudden drop in productivity occurring when going from the latter to the former. Such a discontinuity is well explained by equilibrium theory in terms of operating regimes of humidifier and dehumidifier. Finally, equilibrium theory is used to derive an analytical approximation of the optimal operating curve. Such an analytical form describes well the actual solution, particularly within the boundary conditions of interest, and provides an immediate tool for easily determining the optimal system operation. Keywords: Equilibrium theory, solar-driven desalination, humidification-dehumidification, solar modules, HDH, PVT

1. Introduction Humidification-dehumidification (HDH) systems proved to be a promising solution for small-scale desalination, appropriate for water production in off-grid locations where the water demand does not justify the installation of conventional large-scale systems [1–4]. The HDH process uses air as the carrier fluid to separate clean water from the saline feed, and in its simplest form it consists of three sub-systems: (i) a dehumidifier, where saline water and moist air enter at low and high temperature, respectively; (ii) a humidifier, where saline water and moist air enter at high and low temperature, respectively; (iii) an energy source, which connects the two aforementioned modules. In fact, several configurations have been investigated and ∗ Corresponding

author: [email protected] 2019/07/10

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Hot saline water

Page 2 of 29

PVT modules

3

2 Electricity

Electricity export

Humidifier

Air with high water content 6

Dehumidifier

Packed column

Plate

5

7

Air with low water content

1

Air fan Pure water 4

Water pump Brackish water

Brine

Figure 1: Schematic of a closed-air, open-water, water-heated HDH system with PVT modules as the energy source. Figure adapted from Energy Conversion and Management: X, Vol 1, Gabrielli et al., Combined water desalination and electricity generation through a humidification-dehumidification process integrated with photovoltaic-thermal modules: Design, performance analysis and techno-economic assessment, 100004 [7]. Copyright © 2019 with permission from Elsevier.

reviewed by M¨ uller-Holst et al. [5] and by Narayan et al. [2]. Such studies indicated closed-air, open-water (CAOW), water-heated process as the most energy-efficient solution, due to the high heat capacity of water and the possibility of working with saturated air. A review of the design strategies currently adopted for CAOW processes was recently presented by Elmutasim et al. [6]. A schematic of the HDH system considered in this work, featuring a single-stage CAOW configuration and photovoltaic-thermal solar panels (PVT) as energy source, is shown in Fig. 1. The inlet saline water stream enters the dehumidifier, where it cools down the stream of hot saturated air causing some of the water vapor to condense, thus yielding a stream of clean water. After being pre-heated in the dehumidifier, the saline water enters the solar modules and is heated up to the maximum process temperature. Then, it enters the humidifier, where it heats up the stream of cold saturated air causing some water to evaporate, hence increasing its water vapor content. Finally, after the humidifier the saline water is discharged in the form of concentrated brine. It is worth mentioning that the low productivity of the HDH process, typically lower than 10%, leads to levels of brine concentration, i.e. at the discharge, that allow avoiding expensive brine disposal. Design and optimal operation of the HDH process have been investigated in the past. First, M¨ uller-Holst and Mistry et al. showed that minimizing the entropy generation leads to the minimum energy consumption, and that most of the entropy generation results from the heat and mass transfer in the dehumidifier and humidifier [8–10]. Narayan et al. observed that the entropy generation of combined heat and mass transfer devices is minimized when the heat capacity rates of air and water streams are equal in the dehumidifier [11]. This result leads to an optimal ratio of mass flow rates determined by the operating temperatures of the dehumidifier. Later, Thiel and Lienhard found that the higher relevance of the dehumidifier in the whole process scheme results from the high diffusive mass transfer [12], whereas Chehayeb et al. related it to its higher entropy generation [13]. Based on such findings, several studies presented enhanced layouts of HDH process, where the entropy generation was minimized by exploiting multi-stage configurations to attain the optimal mass flow rate ratio across the entire process, e.g. [13–26]. Nevertheless, it appears that the asymmetry between humidifier and dehumidifier, and particularly the 2

ACS Paragon Plus Environment

Page 3 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

role of the latter in determining the optimal operating conditions of the system, is not well understood, yet. Motivated by this consideration, we provide an alternative perspective with respect to that of references [8–13], which is based on the so-called equilibrium theory, and allows shedding brighter light on the system behavior. Equilibrium theory is applied when dealing with one-dimensional, nonstationary processes that are encountered in various areas of engineering, particularly in separation process technology, such as chromatography, ion-exchange, gas/oil displacement, distillation, and sedimentation but also heat storage, precipitation/dissolution, and traffic flow [27]. Such diverse systems share the same mathematical description, i.e. a system of hyperbolic first-order partial differential equations obtained by considering only convective transport, hence neglecting dispersive effects, and assuming that transport between phases is so fast that thermodynamic equilibrium prevails everywhere in the transfer unit at any time (hence the expression equilibrium theory) [28]. The importance of equilibrium theory stems from four features. First, many different processes of interest can be considered and their analogies emphasized. Second, the key features of the complex dynamic behavior of such systems are explained and predicted accurately. Third, the mathematical tools to tackle this type of equations are complex to develop but simple to apply and more insightful when compared to numerical methods. Finally, equilibrium theory can be used very effectively in the design of complex processes. The underlying idea is that, despite the aforementioned simplifications, equilibrium theory provides accurate predictions and valuable insight into physical processes [27]. Concerning the aforementioned desalination application, we present the equilibrium theory solution valid when assuming a constant saline water flow rate across the humidifier. While this solution assumes the same temperature for the air and water streams within humidifier and dehumidifier, it allows explaining the asymmetry between the two units, which have different features in terms of properties of the temperature and composition fronts that are present in each of them. Within this framework, the relevance of the dehumidifier in determining the optimal operating condition, namely the mass flow rate ratio, results from the functional dependence of the enthalpy of moist air on temperature and from the nature of the dehumidifier, where moist air is cooled down and water is heated up. Furthermore, the derived solution allows for an immediate description of the operation of the overall system, with solar modules connecting the dehumidifier and humidifier, in terms of mass flow rate ratio and solar energy input. The optimal operation curve, namely the mass flow rate ratio that maximizes the productivity for any given ambient conditions, and a robust optimal operation curve are determined exactly. The paper is structured as follows. Section 2 introduces the equilibrium theory equations for a countercurrent heat and mass exchanger dealing with a single-component separation and explains the different types of transitions that can be observed in the solution of the model equations. Such equations are then applied to the dehumidifier and humidifier components in Section 3, whereas the analysis of the overall solar-driven desalination system is presented and discussed in Section 4. Finally, major findings are summarized and conclusions are drawn in Section 5. 2. Equilibrium theory solution for a countercurrent heat and mass exchanger Equilibrium theory solutions have been extensively discussed for a variety of engineering applications [27, 28]. Here we follow the mathematical derivation proposed by Rhee et al. for the analysis of a countercurrent adsorber [28], and we adapt it to describe a countercurrent heat and mass exchanger where a liquid and a gas phase, namely water and air, are contacted. This implies considering a one dimensional convective flow, as well as neglecting dispersive effects and transport resistances. Overall, the following assumptions are made: i) One dimensional convective flow and negligible dispersive effects and transport resistances. ii) Negligible pressure drop and heat losses (i.e. adiabatic operation). iii) Negligible variation of liquid flow rate. iv) Negligible effect of salt content on the properties of saline water.

3

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 29

v) Saturated air across the entire exchanger, i.e the humidity ratio is a known function of temperature, w = wsat (T ), through the vapor pressure of water. More specifically, the humidity ratio is expressed as: w(T ) =

ML pv (T ) MG (p − pv (T ))

(1)

where M is the molar weight and the subscripts ”L” and ”G” indicate water and dry air, respectively; p is the air pressure and pv the vapor pressure of water, which is calculated using Antoine’s equation: log10 (pv ) = C1 −

C2 C3 + T

(2)

where C1 , C2 and C3 are water-specific constants. vi) Constant thermodynamic parameters and volume fractions. Based on such assumptions, the heat and mass exchanger can be described by the energy balance only, that includes the accumulation term and the convective term. More specifically, the liquid and gas streams are characterized in terms of specific mass hold up, mass flow rates and specific enthalpy: liquid gas

specific mass hold up [kg/m3 ]

mass flow rate [kg/s]

specific enthalpy [kJ/kg]

ρL (1 − )ρG

L G

cL (T − Tref ) c∗U (T − Tref )

A schematic of the heat and mass exchanger is shown in Fig. 2. This represents both the humidification and the dehumidification units of the process investigated in this study. The following energy balance can be written: A

∂ ∂ [ρL cL (T − Tref ) + (1 − )ρG c∗U (T − Tref )] + [LcL (T − Tref ) − Gc∗U (T − Tref )] = 0 ∂t ∂z

(3)

where A is the cross sectional area,  the volume fraction occupied by the liquid phase and assumed to be constant, ρ the density, w the air humidity ratio (water content per unit dry air), L and G the mass flow rates of saline water and dry air, respectively, c the specific heat capacity at constant pressure, T the fluid temperature, which is equal for both phases due to the equilibrium assumption, and Tref a reference temperature, which here is chosen equal to 273.15 K; the subscript ”U” indicates moist air, whereas the superscript ”*” refers to a quantity expressed per unit dry air. Here, water flows in the direction of z, whereas air flows in the opposite direction. Note that the mass flow rate of dry air, G, is constant across the separation unit, with the air flow rate variation being caused by the variation of humidity ratio, w. The specific heat capacity of moist air consists of two contributions, namely that of the dry air and that of the water vapor:   h ∗ cU = cG + w(T ) cV + (4) T − Tref a, A liquid in

moist gas out

𝑇a

𝑇A , 𝑓A

b, B

𝑇0+

liquid phase

𝑓0+

moist gas phase

𝐿, 𝜖, 𝜌L , 𝑐L , 𝑇

𝑇1−

𝐺, 1 − 𝜖 , 𝜌G , 𝑐G , 𝑤, 𝑇, 𝑓

𝑓1−

𝑇B

𝑇b , 𝑓b

liquid out

moist gas in

𝑥=1

𝑥=0

Figure 2: Schematic and notation for analysis of heat and mass exchangers. The states just inside the boundaries, denoted by superscripts ”+” and ”-”, are highlighted for both phases.

4

ACS Paragon Plus Environment

Page 5 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

where cV is the specific heat capacity at constant pressure of water vapor and h its latent heat of vaporization. This implies assuming that water is evaporated at Tref and then the water vapor is heated (or cooled) to the gas temperature, T . The phase capacity ratio, ν, and the mass flow rate ratio, µ, are defined as: ν=

(1 − )ρG , ρL

µ=

G L

(5)

with the latter being varied to control the operation of the heat and mass exchanger, as discussed in the following sections. By using such dimensionless parameters, and by recalling the assumptions above, Eq. (3) is recast as: AρL as:

∂ ∂ [cL (T − Tref ) + νc∗U (T − Tref )] + L [cL (T − Tref ) − µc∗U (T − Tref )] = 0 ∂t ∂z

(6)

Furthermore, the dimensionless time coordinate, τ , and the dimensionless space coordinate, x, are defined τ=

tL , AZρL

x=

z Z

(7)

where Z is the length of the exchanger. Moist air is characterized by the function f (T ), which expresses the ratio between the enthalpy of saturated air and the specific heat capacity of water, and has the dimensions of a temperature. It is defined as: f (T ) =

cG (T − Tref ) + w(T ) [cV (T − Tref ) + h] cL

(8)

It plays a key role in determining the solution. Indeed, moist air is characterized by its temperature, T , and by its water content, w. Thus, within the framework of equilibrium theory we characterize it by referring to either T or f (T ), or both. Hence, the following first-order, quasi-linear, homogeneous partial differential equation is obtained: [1 + νf 0 (T )]

∂T ∂T + [1 − µf 0 (T )] =0 ∂τ ∂x

(9)

where f 0 (T ) indicates the first-order derivative of f with respect to T . Note that Eq. (9) presents the same mathematical structure as in distillation, or nonlinear moving-bed adsorption and chromatography, where the function f corresponds to the adsorption isotherm (which is system specific, being dependent on the adsorbent and the adsorbate) and such equilibrium theory model has broad and very effective application [27, 28]. Also, note that the same equation holds for the analysis of a co-current heat and mass exchanger, where the liquid and gas streams have the same direction (the sign ”-” on the second term of the left-hand side is replaced by the sign ”+”). In the case of moist air, f is a universal nonlinear convex function of T with positive first- and second-order derivatives, as shown in Fig. 3-a. A very common and useful application of equilibrium theory is to the so-called Riemann problems. In this case Eq. (9) is solved by considering uniform initial state of the column and constant feed temperatures from both sides. In other words, we solve an initial value problem defined by piecewise constant initial data, where the initial value problem is characterized by three segments and two discontinuities along the initial line, in (0, 0) and (1,0): τ =0 x=0 x=1

0 ≤ x ≤ 1, τ ≥ 0, τ ≥ 0,

T = T0 T = Ta f = fb

f0 = f (T0 ) fa = f (Ta ) Tb = f −1 (fb )

where the subscripts ”0” refer to the initial conditions, ”a” and ”b” to the column ends where there is the liquid inlet (x = 0) and the gas inlet (x = 1), respectively. 5

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 29

Figure 3: a) Function f (T ) describing the equilibrium correlation between the enthalpy of the gas and of the liquid phases. Convex function with positive first- and second-order derivatives. Tref = 273.15 ◦C, salt content of inlet saline water of 20 g/l. b) Illustration of characteristics on the physical plane (x, τ ) for a specific value of temperature T . The characteristic rotates counterclockwise as µ increases, and it is vertical for µ = 1/f 0 (T ).

This type of problems can be solved by using the method of characteristics. For this class of equations, the characteristics are straight lines in the physical plane (x, τ ) that propagate the state on the initial line unaltered [28]. In the case of Eq. (9), the characteristics are lines of constant value of T with slope, σ, in the plane (x, τ ) given by σ(T ) =

dτ 1 + νω = dx 1 − µω

(10)

where ω = f 0 (T ). The slope σ(T ) is positive or negative depending on the value of µ, i.e. of dry air to water mass flow rate ratio, and exhibits a vertical asymptote. The reciprocal of the slope, λ, gives the propagation rate of the state T : λ(T ) =

1 1 − µω = σ(T ) 1 + νω

(11)

which is positive or negative depending on the value of µ, with the vertical asymptote of σ corresponding to λ = 0. This is illustrated in Fig. 3-b, which shows the characteristics on the physical plane (x, τ ) for a given value of T and for different ranges of µ. Such characteristics rotate counterclockwise as µ increases, and they are vertical for µ = 1/ω = 1/f 0 (T ). As µ increases, states slow down in algebraic terms (f 0 (T ) > 0) since: dλ −ω = 0): dλ −(ν + µ)f 00 (T ) = ωR , σL > σR , λL < λR . This situation is illustrated on the left-hand side of Fig. 4. Here, there is a region between the characteristics with slope σ(TL ) (on the left of which there is the constant state TL ) and σ(TR ) (on the right of which there is the constant state TR ) occupied by a so-called simple wave, consisting of a pencil of characteristics emanating from the discontinuity having intermediate slopes and

Figure 4: Characteristics describing the solution on the dimensionless physical plane (x, τ ) and profile of exchanger temperature at x = 1. Left: solution characterized by a simple wave, with T in > T0 . Right: solution characterized by a shock, with T in < T0 . Figure adapted with permission from reference [27]. Copyright © 2013 by Annual Reviews.

7

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 29

propagating intermediate states [29]. Such states increase from TR to TL , as the characteristics fan counterclockwise. The slope of the characteristics varies continuously as described by Eq. (10). ii) TL < TR . Thus, ωL < ωR , σL < σR , λL > λR . This situation is illustrated on the right-hand side of Fig. 4. Here, the triangular region occupied by the simple wave in the previous case would be covered by overlapping characteristics (the left state would catch up with the right state). This leads to a situation in which more than one state would be associated to the same point in the physical plane, in other words, a physically unfeasible solution. This is equivalent to stating that the left state cannot overcome the right state, though the former is faster than the latter. This results in the formation of a shock, i.e. a discontinuity emanating from the initial discontinuity and propagating as a straight line with slope and propagation rate given by Eqs. (10) and (11), with the slope of the tangent to f (T ) being replaced by the chord, i.e. by replacing ω with ω ˜ defined as [28] ω ˜ L,R =

[f ] f (TL ) − f (TR ) = TL − TR [T ]

(15)

where [·] represents the jump of the quantity within brackets across the discontinuity. In the special case where the function f (T ) is linear, then the characteristic slope is independent of the temperature. With reference to our problem, this implies that the characteristics carrying the state TL are parallel to those carrying the state TR . The two corresponding constant states are separated by a straight line of the same slope, which is called a contact discontinuity. 3) Behavior of outgoing streams. There is no reason why the outgoing streams should be in equilibrium with the incoming streams, even if assuming equilibrium inside the exchanger [28]. In fact, the outlet temperatures depend on the contact between the two phases at the boundaries. Assuming no accumulation at the column ends, and using the notation in Fig. 2, one can write the following energy balances: x = 0, x = 1,

Ta − T0+ µ TB = T1− + µ(fb − f1− )

fA = f0+ +

(16) (17)

which determine the outlet temperature of the gas phase, TA , (through the function fA = f (TA )) and the outlet temperature of the liquid phase, TB . Here, the subscripts ”+” and ”-” represent the states just inside the boundaries at x = 0 and x = 1, respectively; the lower case subscripts, ”a” and ”b” refer to the inlet conditions, whereas the upper case subscripts, ”A” and ”B” refer to the outlet conditions at the corresponding column ends. Therefore, Ta and fb are known from the formulation of the initial value problem, whereas the states T and f with subscripts ”0+” and ”1−” are obtained from the solution of the PDE. These boundary conditions are then used to determine the outlet states once the states inside the exchanger are known. In this case, temperature discontinuities are observed at both boundaries. This type of discontinuities is different than either a shock or a contact discontinuity, and is called boundary discontinuity [28]. One can note that such discontinuities imply an energy loss, and the ideal solution from a thermodynamic perspective is such that no boundary discontinuity occurs. For the sake of simplicity, we now illustrate the solution of the heat and mass exchanger by considering a simplified Riemann problem with Tb = T0 (initial state and gas feed state - entering from the right - are the same), hence with only one discontinuity in (0, 0). We want to study the structure of the solution when varying µ, i.e. the manipulated operating parameter, and we need to consider the same two cases introduced above, namely simple waves and shocks as possible forms of the solution. i) Ta > T0 . The solution is described by a simple wave centered at the origin, with three possible situations occurring based on the value of µ: 8

ACS Paragon Plus Environment

Page 9 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

• For low values of µ the simple wave consists of characteristics propagating forward, and the solution on the physical plane (x, τ ) looks like that shown in Fig. 5-a. The structure of the solution remains unchanged as long as the characteristic associated to Ta is directed forward, thus for 0 ≤ µ < µ∗ =

1 f 0 (Ta )

(18)

which is called the lower range of µ. Within this range T0+ = Ta , hence fA = f0+ = fa , i.e. there is no discontinuity at x = 0 and equilibrium is maintained across the boundary. On the other hand, the nature of the boundary at x = 1 is different and varies with time. By calling τ0 = σ(T0 ) and τa = σ(Ta ), the outlet liquid temperature is given by 0 ≤ τ < τ0 ,

TB = T0

τ0 ≤ τ < τa ,

TB = T + µ(f0 − f ),

τa ≤ τ ,

TB = Ta + µ(f0 − fa )

Ta ≥ T ≥ T0

For τ < τ0 , the temperature at x = 1 is the initial one as the front of the simple wave has not reached there yet. For τ0 ≤ τ < τa , T increases from T0 to Ta ; equilibrium is not maintained across the boundary and TB is determined through Eq. (17). For τ ≥ τa , no further change is observed, and the system reaches the following steady state: 0 < x < 1, T0+ = T1− = Ta x = 0,

fA = fa

x = 1,

TB = Ta + µ (f0 − fa )

(19)

with the state Ta occupying the whole column. A better feeling of the variation of TB is obtained by referring to the function f (T ), as illustrated in Fig. 5-b, where the boundary discontinuity at x = 1 is represented by the line segment of slope 1/µ. It is clear that TB decreases with time; in particular, it follows the profile shown in Fig. 4 for the simple wave.

Figure 5: a) Characteristics describing the solution on the dimensionless physical plane (x, τ ) and (b) profile of exchanger temperature just inside the boundary at x = 1 for 0 ≤ µ < µ∗ . Solution described by a simple wave centered in the origin with all characteristics directed forward (blue solid lines), with a boundary discontinuity at x = 1 (black dashed line); τ0 and τa denote the dimensionless time needed for the simple wave to reach x = 1 and for the system to reach the steady state, respectively. Operating parameters: T0 = 303 K, Ta = 318 K, ν = 0.011, µ = 0.3, µ∗ = 0.37.

9

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 29

• The structure of the solution changes when µ is increased to fall within the range 1 1 = µ∗ ≤ µ < µ∗ = 0 f 0 (Ta ) f (T0 )

(20)

which is called the transition range of µ (see Fig. 6). Within this range, characteristics have partly positive slope (as shown in Fig. 6-a) and partly negative slope (not seen in Fig. 6-a as they propagate backward from the vertical axis x = 0); hence, there is one state, T∞ such that Ta ≤ T∞ < T0 , whose characteristic is vertical: f 0 (T∞ ) = µ1 ; this is the only state to remain present within the exchanger at t → ∞. Therefore, the steady state is only reached at infinite time, where: 0 < x < 1, T0+ = T1− = T∞ Ta − T∞ µ

x = 0,

fA = f∞ +

x = 1,

TB = T∞ + µ (f0 − f∞ )

(21)

• Finally, the structure of the solution changes again when further increasing µ, to reach the upper range of µ, µ ≥ µ∗ =

1 f 0 (T0 )

(22)

In this case, the characteristics for all temperatures between the initial and the inlet values are directed backward. Therefore, the state at x = 0 cannot penetrate the boundary, where a boundary discontinuity occurs, while the state at x = 1 is propagated into the exchanger, where equilibrium

Figure 6: a) Characteristics describing the solution on the dimensionless physical plane (x, τ ) and (b) profile of exchanger temperature just inside the boundary at x = 1 for µ∗ ≤ µ < µ∗ . Solution described by a simple wave centered in the origin with characteristics directed both forward and backward (blue solid lines), with boundary discontinuities at x = 0 and x = 1 (black dashed line); Operating parameters: T0 = 303 K, Ta = 318 K, ν = 0.011, µ = 0.6, µ∗ = 0.37, µ∗ = 0.82.

10

ACS Paragon Plus Environment

Page 11 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

is maintained across the boundary. The system attains the following steady state: 0 < x < 1, T0+ = T1− = T0 x = 0,

fA = f0 +

x = 1,

TB = T0

Ta − T0 µ

(23)

ii) Ta < T0 . A shock propagates from the origin instead of the simple wave discussed above. The shock maintains a constant strength of (Ta − T0 ) and propagates in the constant direction given by σ ˜a,0 =

1 + νω ˜ 1 − µ˜ ω

(24)

where ω ˜=

fa − f0 Ta − T0

(25)

In this case, the transition range of µ shrinks to a critical value µc , corresponding to µc = µ∗ = µ∗ =

1 ω ˜

(26)

for which the following steady state is reached: 0 < x < 1, T0+ = Ta , f1− = f0 x = 0,

fA = fa

x = 1,

TB = T0

(27)

Such condition is optimal from the perspective of energy exchange, as both outgoing streams are at equilibrium with the corresponding incoming streams (with T0 = Tb ). In general, the initial state of the mass and heat exchanger is different than the two inlet states, liquid from x = 0 and gas from x = 1. Thus, two waves must be considered, one emanating from the origin and the other from the point (1, 0) on the plane (x, τ ). If at least one of the two is a simple wave, the analysis is in principle similar to what discussed above. However, the two simple waves may meet inside the exchanger and interact; the description of such interaction is possible but complex (see [28]). It is however outside the scope of this work, as these interactions do not affect the steady state behavior of the HDH process. On the contrary, if a shock occurs at both sides, and if they meet inside the exchanger and interact with each other, such interaction can be analyzed in a rather simple way. This situation is discussed in the following for the humidification-dehumidification process. Before diving into this, we clarify why the full problem, and not the steady state version of it only, must be solved also when investigating the steady state of the system. Because of the homogeneous nature of Eq. (6), any solution that is constant along the space coordinate would fulfill such equation. One can readily see that, for the type of Riemann problem that we are interested to solve in practice, there is no physical reason to select one or the other temperature level, as long as we take one within the inlet temperatures of the gas and liquid phases. The only steady states that are consistent with the original dynamic nature of the problem can be selected by considering the dynamic nature of the original problem. This allows explaining the asymmetry between the humidification and the dehumidification processes, which is based on the different propagation of the temperature and composition fronts within the two exchangers, and cannot be captured by solving the steady-state model given by Eqs. (16) and (17) only. 11

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 29

3. Equilibrium theory solution for a humidification-dehumidification process In this section, the equilibrium theory model is applied to study the humidifier and the dehumidifier components of the HDH process, first separately, and then in the proper connected configuration of Fig. 1. Before this, the special case of a countercurrent heat exchanger with no mass transfer is investigated. While the behavior of the exchangers is investigated on the physical plane (x, τ ), the analysis focuses in particular on the steady state of the system, which depends only on the value of µ and is independent of the initial state. 3.1. Countercurrent heat exchanger The special, simple case of a countercurrent heat exchanger is investigated to illustrate the comparison between the equilibrium theory and the results obtained by Narayan et al. in Section 2 of reference [11]. ref ) . Eq. (9) is written as Here, dry air is considered, i.e. w = 0 and f = cG (Tc−T L (1 + νω)

∂T ∂T + (1 − µω) =0 ∂τ ∂x

(28)

where ω = cG /cL is a constant. This implies that all characteristics have the same slope, which is independent of temperature and only function of µ: σ=

1 + νω 1 − µω

(29)

Therefore, the solution consists of parallel characteristics that carry the initial and the inlet states when emanating from the x- and τ -axis, respectively, and are separated by a contact discontinuity. The solution is shown in Fig. 7 for the different ranges of µ, separated by the critical value µc = 1/ω. Within the lower range of µ, 0 ≤ µ < µc , all characteristics are directed forward and the heat exchanger is flooded by the state of the liquid phase, reaching the temperature of the incoming water at steady state. Within the upper range of µ, µ > µc , all characteristics are directed backward and at steady state the heat exchanger is at the temperature of the incoming gas. When µ = µc , all characteristics are vertical, being σ infinitely large. Thus, we observe vertical contact discontinuities at both boundaries, i.e. the inlet states of both the liquid

Figure 7: Equilibrium theory solution on the plane (x,τ ) for a countercurrent heat exchanger within different ranges of µ. Contact discontinuity indicated by the green line. Operating parameters: ν = 0.011, ω = 0.24, (a) µ = 2, 0 ≤ µ < 1/ω, (b) µ = 4.16, µ = 1/ω, (c) µ = 6, µ > 1/ω.

12

ACS Paragon Plus Environment

Page 13 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Table 1: Steady state of countercurrent heat exchanger for different ranges of µ.

0 < x < 1 (exchanger)

x = 0 (inlet liquid)

x = 1 (inlet gas)

µ < µc

T = Ta

fA = fa

TB = Ta + µ (fb − fa )

µ = µc

T = Ta |Tb

fA = fa

µ > µc

T = Tb

fA = fb +

TB = Tb 1 µ

(Ta − Tb )

TB = Tb

and the gas phase remain stationary at x = 0 and x = 1, respectively. Such discontinuities are present from the beginning and the steady state is reached immediately, with equilibrium being maintained across both boundaries. The steady states attained by the system for the different ranges of µ are reported in Table 1, where the notation Ta |Tb indicates that both states Ta and Tb are present inside the exchanger, i.e. T0+ = Ta and T1− = Tb . Based on such temperature profiles, the entropy generation of the heat exchanger, S, is calculated as     TA TB + GcG log (30) S˙ = LcL log Ta Tb To compare the findings of equilibrium theory to those obtained by Narayan et al. [11], we compute the dimensionless entropy generation, s, ˙ by dividing S by the minimum heat capacity rate of liquid and gas phase, i.e. GcG for µ < µc and LcL for µ > µc . Replacing the corresponding expressions in Table 1, s˙ is obtained as a function of µ as follows:     1 µω(Tb − Ta ) Ta µ < µc , s˙ = log 1 + + log ωµ Ta Tb     Tb Ta µ = µc , s˙ = log + log =0 (31) Ta Tb

.

Figure 8: Dimensionless entropy generation, normalized over the minimum heat capacity rate of liquid and gas phase, as function of µ. Ta = 2Tb .

13

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

s˙ = log

µ > µc ,



Tb Ta



Page 14 of 29

  Ta − Tb + µω log 1 + µω(Tb )

The results are illustrated in Fig. 8, which shows that the entropy generation of the heat exchanger is minimized, hence the efficiency of the heat exchanger is maximized, when µ = µc . This condition is equivalent to stating that the heat capacity ratio must be equal to one, i.e. HCR = 1, as determined by Narayan et al. [11]: µω = HCR and µc = 1/ω. In this case, both outgoing streams are at the same temperature of the corresponding incoming streams, and no entropy is generated in the process. While leading to an entropy production equal to zero, the equilibrium theory solution is able to determine (i) the condition of minimum entropy generation, (ii) the qualitative behavior of the generated entropy for values of µ different from µc . The quantitative comparison between our results and those obtained by Narayan et al. [11] is carried out in Section S1 of the Supporting Information, where we introduce the concept of heat exchanger effectiveness, which allows accounting for heat transfer limitations starting from the equilibrium theory solution. 3.2. Dehumidifier The dehumidification unit of the desalination process is now investigated, with the equilibrium theory solution being studied for the case of the incoming liquid being always colder than the incoming gas, i.e. Ta < Tb . The types of transitions connecting the contact states Ta , T0 and Tb in the physical plane depend on the initial temperature level, T0 , relative to Ta and Tb . The ultimate steady-state depends only on the inlet states Ta and Tb . More specifically, based on the values of Ta , T0 and Tb the following transitions, associated to the initial discontinuities in (0,0) and in (1,0), appear in the (x, τ ) plane:

Ta < Tb < T0 Ta < T0 < Tb T0 < Ta < Tb

Inlet water (0,0)

Inlet moist air (1,0)

shock shock simple wave

simple wave shock shock

We illustrate the behavior of the dehumidifier by referring for simplicity to the case of two shocks, where Ta < T0 < Tb . Here, the form of the solution depends on the value of µ and on the following characteristic slopes: σ ˜a,0 =

1 + νω ˜ a,0 , 1 − µ˜ ωa,0

σ ˜0,b =

1 + νω ˜ 0,b , 1 − µ˜ ω0,b

σ ˜a,b =

1 + νω ˜ a,b 1 − µ˜ ωa,b

(32)

fa − fb Ta − Tb

(33)

with, ω ˜ a,0 =

fa − f0 , Ta − T0

ω ˜ 0,b =

f0 − fb , T0 − Tb

ω ˜ a,b =

More specifically, three possible steady states may occur: (i) µ < µc , exchanger occupied by the liquid phase, (ii) µ = µc , both phases present inside the exchanger, (iii) µ > µc , exchanger occupied by the gas phase [28], where the critical value µc is defined as µc =

1 Ta − Tb = ω ˜ a,b fa − fb

(34)

Such steady states are described below. i) µ < µc , lower range of µ. The shock from the origin propagates forward since σ ˜a,0 > 0. Different situations are observed concerning the system transient: • For small values of µ, µ < 1/˜ ω0,b , the shock emanating from the gas side does not come into the exchanger at all, as the slope σ ˜0,b is positive. 14

ACS Paragon Plus Environment

Page 15 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

• As µ increases above this value, the slope σ ˜b,0 becomes negative and the shock propagates backward into the exchanger. The two shocks then meet, combine, and continue as a single shock. This phenomenon is called the confluence of shocks [28, 30]. The combined shock of strength (Ta − Tb ) propagates in the direction σ ˜a,b > 0 and eventually reaches the boundary at x = 1. Thereafter, the system reaches the steady state given by Eq. (19), with T0+ = T1− = Ta and a discontinuity at x = 1. This situation is depicted in Fig. 9-a. ii) µ = µc , critical value of µ. The two shocks meet and combine as before, but in this case the combined shock becomes stationary because σ ˜a,b → ∞. Thanks to the presence of the standing shock wave, there are no discontinuities at either boundary, so that equilibrium is maintained at both sides of the exchanger. This situation is illustrated in Fig. 9-b and -d on the physical plane and in terms of f (T ),

Figure 9: Equilibrium theory solution for a countercurrent air dehumidifier for different ranges of µ. Shocks indicated by the red lines. Operating parameters: Ta = 293 K, Tb = 318 K, T0 = 313 K, ν = 0.011, µc = 0.67, (a) µ = 0.5, 0 ≤ µ < µc , (b) µ = 0.67, µ = µc , (c) µ = 0.8, µ > µc . d) Profile of exchanger temperature just inside the boundary at x = 1 for µ = µc .

15

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 29

respectively. The attained steady state is given by: 0 < x < 1, T0+ = Ta , T1− = Tb x = 0,

fA = fa

x = 1,

TB = Tb

(35)

This is the best result that can be achieved in terms of heat and mass exchange, as both outgoing streams reach equilibrium with the corresponding incoming streams and no energy loss occurs. iii) µ > µc , upper range of µ. The combined shock (or only the shock emanating from the gas side for µ>ω ˜ a,0 ) has a negative slope and thus propagates backward to reach the boundary at x = 0, as shown in Fig. 9-c. Here, T0+ = T1− = Tb for 0 < x ≤ 1, with a boundary discontinuity at x = 0, which is the same steady state given by Eq. (23). The steady states of the dehumidifier are summarized in Table 2 together with those of the humidifier. As for Table 1, the condition for which T0+ = Ta and T1− = Tb is indicated as Ta |Tb . 3.3. Humidifier The humidification unit of the desalination process is now investigated, with the equilibrium theory solution being studied for the case of the incoming liquid always warmer than the incoming gas, i.e. Ta > Tb . As for the dehumidifier, the types of transitions occurring during the transient depend on the initial exchanger temperature, T0 , whereas the ultimate, steady state transitions depend only on the inlet states Ta and Tb . More specifically, based on the values of Ta , T0 and Tb the following transitions occur on the (x, τ ) plane:

Ta > Tb > T0 Ta > T0 > Tb T0 > Ta > Tb

Inlet water (0,0)

Inlet moist air (1,0)

simple wave simple wave shock

shock simple wave simple wave

Contrary to the dehumidifier, the case of two shocks occurring at both sides of the exchanger is not observed, and at least one simple wave is always present inside the exchanger. In this case, the analysis of the interaction is more complex and we discuss it in Section S2 of the Supporting Information with reference to the case of T0 > Ta > Tb , i.e. with a shock propagating from the origin and a simple wave propagating from (1, 0). In general, the form of the steady state depends on the lower and upper limit values µ∗ and µ∗ (see Section 2): µ∗ =

1 , fa0

µ∗ =

1 , fb0

(36)

Table 2: Steady state of countercurrent humidifier and dehumidifier for different ranges of µ.

Dehumidifier (Ta < Tb )

T (0 < x < 1)

fA (x = 0)

TB (x = 1)

µ < µc

Ta

fa

Ta + µ (fb − fa )

µ = µc

Ta |Tb

fa

Tb

fb +

Ta

fa

Ta + µ (fb − fa )

T∞

f∞ + µ1 (Ta − T∞ ) fb + µ1 (Ta − Tb )

T∞ + µ (fb − T∞ )

µ > µc Humidifier (Ta > Tb )

µ < µ∗ µ∗ ≤ µ < µ ∗

µ≥µ



Tb

Tb 1 µ

(Ta − Tb )

16

ACS Paragon Plus Environment

Tb

Tb

Page 17 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

More specifically, three possible steady states may occur: (i) µ < µ∗ , exchanger occupied by the state of the liquid phase, (ii) µ∗ ≤ µ < µ∗ , exchanger occupied by an intermediate state T∞ , with Tb < T∞ < T0 and f 0 (T∞ ) = 1/µ, which is attained at infinite time (see Eq. (21)), (iii) µ > µ∗ , exchanger occupied by the state of the gas phase [28]. 3.4. Steady state of humidifier and dehumidifier The steady states of humidifier and dehumidifier for the specific inlet temperatures considered above are also shown in Fig. 10 in terms of outlet temperatures, dimensionless entropy generation, and evaporated/condensed clean water per unit dry air, W . The evaporated water in the humidifier is calculated as

Figure 10: Steady state of humidifier/dehumidifier (left/right) as function of µ. Top: outlet temperatures. Middle: dimensionless entropy generation, normalized over the minimum heat capacity rate. Bottom: evaporated/condensed water per unit dry air. Operating parameters: Ta = 318/293 K, Tb = 293/318 K, ν = 0.011, µ∗ = 0.37, µ∗ = 1.4, µc = 0.67.

17

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 29

the difference in humidity ratio between the outlet and inlet moist air streams: W = w(TA ) − w(Tb )

(37)

while the condensed water in the dehumidifier is calculated as the difference in humidity ratio between the inlet and the outlet moist air streams: W = w(Tb ) − w(TA )

(38)

The outlet temperatures (top of Fig. 10) follow the behavior discussed above. Concerning the humidifier, a smooth behavior is observed due to the simple wave interaction between the inlet states of the two phases. For µ < µ∗ the outgoing air is at equilibrium with the incoming water, whereas for µ ≥ µ∗ the outgoing water is at equilibrium with the incoming air. However, there is no value of µ that allows attaining equilibrium at both sides at the same time. Note that the outlet temperatures reported in the interval [µ∗ , µ∗ ] are in fact reached only asymptotically at infinite time. On the contrary, a sharp transition is observed for µ = µc in the dehumidifier, due to the shock interaction between the inlet states of the two phases. In this condition, equilibrium is attained at both sides at the same time. As mentioned above, this represents the optimal condition in terms of heat and mass exchange, and leads to a zero entropy generation when calculating it through Eq. (30). These findings show how the asymmetry between humidifier and dehumidifier does not result from mass transport limitations, but from the convex functional dependence of the enthalpy of moist air on temperature (see Fig. 3), which in turn depends on the functional dependence on temperature of the water vapor pressure. On the one hand, humidifying air becomes easier when increasing the temperature, i.e. smaller temperature variations lead to greater humidity variations, leading to the simple wave interaction and the corresponding smooth transition; on the other hand, dehumidifying air becomes more difficult when decreasing the temperature, i.e. greater temperature variations lead to smaller humidity variations, leading to the shock interaction and the corresponding sharp transition. As a consequence of such transitions, the profile of entropy generation (middle of Fig. 10) is also sharper for the dehumidifier than for the humidifier (see Eq. (30)). This implies that moving away from the optimal condition for the dehumidifier, e.g. working at the optimal condition of the humidifier, would result in a greater variation of entropy generation. Similar considerations apply to the profile of evaporated/condensed water (bottom of Fig. 10). For both units, a plateau at the maximum value is achieved within the lower range of µ, where the outgoing air is at equilibrium with the incoming water. Note that for the two units considered separately, the temperature of the incoming air remains constant, leading to the plateau observed in the figure. However, when investigating the overall HDH system the temperatures of the incoming streams vary based on the interaction between the two units. The coupled system is analyzed in the next section. 4. Analysis of the solar-driven HDH desalination system With the humidification and dehumidification units being characterized in terms of properties of temperature and compositions fronts, we now study the steady-state behavior of the overall desalination process to determine its optimal operating conditions. Here, the humidifier and dehumidifier are coupled through solar modules that heat up the water stream after it exits the dehumidifier and before it enters the humidifier. Hereafter, we refer to the notation in Fig. 11, which shows a scheme of the process where the water and air paths are indicated as blue and red lines, respectively. The saline water enters the dehumidifier at T1 , is preheated within the dehumidifier up to T2 , is heated by the solar panels up to T3 , and is finally cooled down to T4 within the humidifier. On the other side of the exchangers, the air circulates in a closed loop where it is cooled down from T6 to T5 in the dehumidifier and heated up from T5 to T6 in the humidifier. Clean desalinated water is recovered from the dehumidifier. The temperatures inside humidifier and dehumidifier are TH and TD , respectively. The energy input coming from the solar modules is represented by the 18

ACS Paragon Plus Environment

Page 19 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

PVT 𝜇c 𝑇1

𝛼

𝑇2

Dehumidifier 𝑇D

𝜇∗

𝜇∗ 𝑇3

𝑇4

Humidifier 𝑇H

𝑇6 , 𝑓6 𝑇5 , 𝑓5

𝑇5 , 𝑓5

Figure 11: Schematic of the solar-driven HDH process, where the humidification and dehumidification units are coupled by solar modules. The water and air paths are indicated by the blue and red lines, respectively.

lumped parameter α, which is defined below and has the units of a temperature. With this notation, in the dehumidifier the inlet liquid is at T1 and the inlet gas at T6 , with T1 < T6 , hence the critical value of µ is µc =

T1 − T6 f1 − f6

(39)

In the humidifier, the inlet liquid is at T3 and the inlet gas at T5 , with T3 > T5 , hence the two limit values of µ are µ∗ =

1 , f 0 (T3 )

µ∗ =

1 f 0 (T5 )

(40)

The solar modules are accounted for through the energy balance written for a control volume containing the modules themselves, i.e. LcL (T3 − T2 ) = ηIA

(41)

where η is the thermal efficiency of the solar modules, I the solar irradiance and A the installed solar area. In general, the thermal efficiency depends on the ambient conditions and on the temperature of the fluid going through the modules [7]. Here we assume a constant value for η, as this results in an acceptable loss of accuracy (as shown in Section S3 of the Supporting Information) and allows representing the solar modules through a single parameter α, with α = (ηIA)/(LcL ) ≥ 0, having temperature units; α includes both design quantities, such as the specific solar area A/L, and operation quantities, such as the solar irradiance. The overall system is then described by the following system of equations, which includes the energy balance of the solar modules, and the steady state equilibrium theory solutions for the humidifier and the dehumidifier: T3 = T2 + α   0 ≤ µ < µc T1 , TD = T1 |T6 , µ = µc   T6 , µ > µc    T3 ,   0 ≤ µ < µ∗ 0 −1 1 ∗ TH = (f ) µ , µ∗ ≤ µ < µ   T , µ ≥ µ∗

(42)

T2 = TD + µ (f6 − fD )

(45)

T4 = TH + µ (f5 − fH )

(46)

T1 = TD + µ (f5 − fD )

(47)

(43)

(44)

5

19

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

T3 = TH + µ (f6 − fH )

Page 20 of 29

(48)

where Eq. (42) is the energy balance of the solar modules obtained by rearranging Eq. (41); Eqs. (43) to (48) describe the humidification-dehumidification process: Eqs. (43) and (44) describe the humidifier and the dehumidifier units, respectively, Eqs. (45) and (46) determine the outlet liquid temperatures, and Eqs. (47) and (48) determine the outlet gas temperatures. Eq. (44), for µ∗ ≤ µ < µ∗ , means that TH is such that f 0 (TH ) = 1/µ. Note that µ is the same for the two units because of their operational arrangement. The system of Eqs. (42) to (48) consists of 7 equations involving 10 physical quantities, namely T1 , T2 , T3 , T4 , T5 , T6 , TH , TD , µ and α; it has therefore three degrees of freedom. This means that the state of the system can be completely determined by fixing the values of three of these 10 variables, and the remaining 7 unknowns can be calculated by solving the 7 equations. Here, we fix (i) the inlet temperature of the saline water, T1 , which is often defined by the available source of saline water; (ii) α, which specifies the variable ambient conditions under which the desalination system operates; (iii) µ, which is the control variable defining the system operation. The scope of the analysis is to determine the optimal value of the operating variable µ for any set of external conditions given by α and T1 . However, for the coupled system the limit and the critical values of µ, given by Eqs. (39) and (40) depend on the value of µ itself, and must be calculated while determining the solution. Therefore, we first define the relative magnitude of µ∗ , µ∗ and µc and then we determine the solution within the different ranges of values of µ. Before doing this, we note that independently of the values of µ and α Eqs. (42) and (45) to (48) can be rearranged to obtain: TD − TH + α = µ (fD − fH )

(49)

T4 = T1 + α

(50)

where Eq. (49) correlates the states of the humidifier and of the dehumidifier, and Eq. (50) describes the energy balance across a control volume containing the entire system. The latter says that at steady state and assuming that the clean water produced is negligible, the energy input from the solar modules can only increase the enthalpy of the water entering the system at T1 and leaving it at T4 . The relative magnitude of the threshold values of µ is discussed as follows. • µ∗ = 1/f30 defines the lower range of µ, as T3 is the highest temperature of the process, i.e. µ∗ ≤ µc , µ∗ ≤ µ∗ . • The relative magnitude of µ∗ and µc depends on the values of T1 , T5 and T6 . Let us consider the case where µ∗ ≤ µ ≤ µc . Eqs. (43) and (44) state that TD = T1 ,

TH = T5

At the same time, Eqs. (46) and (47) lead to T4 = T5 ,

T5 = T1 =⇒ T4 = T1

whereas Eqs. (45) and (48) lead to T2 = T3 . This is possible only for the trivial solution where α = 0, no energy is provided, and the entire system is at equilibrium at T = T1 . Therefore, for a generic α ≥ 0, µc must be smaller than µ∗ . • µ∗ = µ∗ only if T3 = T5 , which is again possible for the trivial solution only, with no energy provided and the entire system at equilibrium. This implies that µ∗ must be greater than µ∗ , and leads to the following chain of inequalities: µ∗ ≤ µc < µ∗

(51)

After establishing those inequalities, the system of Eqs. (42) to (48) can be rearranged to determine all process temperatures as function of T1 , α, and µ, within the four different ranges of µ, as reported in Table 3. Going from the highest to the lowest value of µ: 20

ACS Paragon Plus Environment

Page 21 of 29

A) µ∗ ≤ µ, upper range of µ. The temperatures inside the dehumidifier and the humidifier are given by Eqs. (43) and (44), respectively: TD = T6 , TH = T5 Thus, T2 = T6 from Eq. (45) and T4 = T5 from Eq. (46). With T1 known, the other temperatures are found by using Eqs. (42), (47) and (48): T4 = T5 = TH = T1 + α T6 − T5 + α = µ(f6 − f5 ) =⇒ T6 , T2 , TD T3 = T2 + α This means that the limit value defining the upper range of µ is determined by T1 and α only, i.e. µ∗ = 1/f 0 (T1 + α). B) µc < µ < µ∗ , upper transition range of µ (approaching µc from above). The temperatures inside the dehumidifier and the humidifier are given by: TD = T6 ,

TH : fH0 =

1 µ

Thus, T2 = T6 from Eq. (45). With T1 known, and TH calculated from the value of µ, the other temperatures are given by: T4 = T1 + α T6 − TH + α = µ(f6 − fH ) =⇒ T6 , T2 T1 = T6 + µ(f5 − f6 ) =⇒ T5 T3 = T2 + α C) µ = µc , critical value of µ. The temperatures inside the dehumidifier and the humidifier are given by: TD = T1 |T6 ,

TH : fH0 =

1 µ

Table 3: Expressions used to determine all process temperatures within the different ranges of µ, and expressions of limit and critical values of µ as function of T1 , α and µ.

range of µ

process variables

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

µ < µ∗

µ = µ∗ = µc

µc < µ < µ∗

µ∗ ≤ µ

TH

T3

fH0 = 1/µ

fH0 = 1/µ

T5

TD

T1

T1 |T6

T6

T6

T2

T3 − α

T6

T6

T6

T3

T1 − T3 + α = µ(f1 − f3 )

TH

T2 + α

T2 + α

T4

T1 + α

T1 + α

T1 + α

T1 + α

T5

T1

T1

f5 = f6 + (T1 − T6 )/µ

T1 + α

T6

T3

µ = (T1 − T6 )/(f1 − f6 )

T6 − TH + α = µ(f6 − fH )

T6 − T5 + α = µ(f6 − f5 )

µ∗

T1 − TH + α = 1/fH0 (f1 − fH ) =⇒ TH =⇒ µ = 1/fH0

µc

T1 − TH + α = 1/fH0 (f1 − fH ) =⇒ TH =⇒ µ = 1/fH0

µ∗

1/f 0 (T1 + α)

21

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 29

Thus, T5 = T1 from Eq. (47) and T2 = T6 from Eq. (45). With T1 known, and TH calculated from the value of µ, the other temperatures are given by: T4 = T1 + α T6 − TH + α = µ(f6 − fH ) =⇒ T6 , T2 T3 = T2 + α D) µ∗ ≤ µ < µc , lower transition range of µ (approaching µc from below). The temperatures inside the dehumidifier and the humidifier are given by: TD = T1 ,

TH : fH0 =

1 µ

Thus, T5 = T1 from Eq. (47). With T1 known, and TH calculated from the value of µ, the other temperatures are given by: T4 = T1 + α T2 = T3 − α = T1 + µ(f6 − f1 ) T3 = TH + µ(f6 − fH )

) =⇒ T2 , T3 , T6

Here Eqs. (42), (45) and (48) must be solved together to determine T2 , T3 and T6 . In general, this cannot be solved since in this case Eq. (49) establishes a relationship between T1 and TH , as discussed more in detail below. E) µ < µ∗ , lower range of µ. The temperatures inside the dehumidifier and the humidifier are given by: TD = T1 ,

TH = T3

Thus, T5 = T1 and T6 = T3 from Eqs. (47) and (48), respectively. With T1 known, the other temperatures are found by: T4 = T1 + α T1 − T3 + α = µ(f1 − f3 ) =⇒ T3 , T6 , TH T2 = T3 − α For both ranges (C) and (D), where TD = T1 and µ = 1/fH0 , Eq. (49) can be written as T1 − TH + α = fH0 (f1 − fH )

(52)

and can be used to calculate TH independently of µ. This means that for each pair (T1 ,α) there is one and only one possible value of TH , hence one possible value of µ only, within the range µ∗ ≤ µ ≤ µc . In other words, such a range (i.e. the lower part of the transition range) shrinks to a single value of µ, i.e. µc = µ∗ = 1/fH0 = 1/f30 = (T1 − T6 )/(f1 − f6 ), which is determined by T1 and α only. One can also see this by looking at the structure of the solution for the µ-ranges in cases (C) and (D): the value of TH can be determined either from the value of µ or from Eq. (49), in fact removing one degree of freedom. On the one hand, this is not surprising when looking at the bottom of Fig. 8: for a coupled system, the value of µ defining the maximum evaporated/condensed water must coincide for the two exchangers. On the other hand, the shrinking of range (D) to the critical value results in a discontinuity at µc , with the system going directly from (B) to (E) when decreasing the value of µ, and vice versa. More specifically, this behavior of the transition range results in sharp profiles of the process temperatures as function of µ and α. This is illustrated in Fig. 12-a, which shows T3 as a function of µ for several values of α. One can see that for each value of α, and for fixed T1 , the corresponding T3 variation as a function 22

ACS Paragon Plus Environment

Page 23 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 12: Illustration of sharp transition range of µ due to µ∗ = µc . (a) Behavior of maximum process temperature, T3 , as a function of µ for several values of α and forfixed T1 . (b) Behavior of solar module-coupled temperatures, T2 and T3 , as function of α for a fixed value of µ = µc α0 , T1 . Temperature discontinuities highlighted. T1 = 273 K, α0 = 4.38, 0 ≤ α ≤ 10.

of µ exhibits a maximum at µ = µc (α): (i) increasing µ above this optimal value (i.e. through the upper transition range and the upper range of µ) leads to a gentle decline in T3 ; (ii) decreasing µ below this optimal value leads to a discontinuity in T3 , with µ entering directly into the lower range. The magnitude of such a jump in temperature increases with increasing α. Similar trends are observed for all other process temperatures. This is illustrated from a different perspective in Fig. 12-b, which shows  the temperatures across the solar modules, T2 and T3 , as function of α for a fixed value of µ = µc α0 , T1 . At α = α0 , which corresponds to the critical value of µ, a discontinuity is observed for both T2 and T3 . The behavior of the temperature profiles impacts directly the productivity, π, of the desalination process, which is defined as the amount of clean water produced per unit saline water treated: π = µ [w(T6 ) − w(T5 )]

(53)

The behavior of the productivity is illustrated in Fig. 13, which shows for a fixed value of T1 the contour lines of π on the plane (µ, α), which is the plane defining the system operation. The contour lines are reported in gray, while the optimal operation line, i.e. µ = µc (α, T1 ), which defines the maximum system productivity, is reported in red. While a smooth decrease in productivity is observed away from the optimal line (both above and below it), a discontinuity in productivity is observed when crossing the optimal line. Ideally, for any ambient condition, µ should vary by moving along the optimal operation line. However, this implies the risk of significant productivity drops for little reductions of µ and α. Therefore, a robust optimal operation would result in a trajectory laying on the right of the optimal line in the (µ, α) plane. In this case, although a lower productivity is achieved with respect to the ideal values of µ, high values of productivity are ensured also for variations of µ and α. An arbitrary robust optimal operation line is illustrated in green, which is built to stay above the optimal operation line in case of a reduction of α of 50%, i.e. µ is selected such that αrob (µ) = 1.5αopt (µ). Finally, we investigate the possibility of determining an analytical expression for the optimal operation line. Although we are not able to find an analytical solution when using the actual f (T ) given by Eq. (8), which is an exponential function of T through the humidity ratio, we can do so by using an empirical approximation of f (T ) given by a convenient functional form, fˆ(T ). Letting θ = T − Tref , and noting that 23

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 29

Figure 13: Operating region of the desalination process and indication of optimal operating conditions. Productivity contour lines are reported in gray, the optimal operation curve is reported in red and the robust optimal operation curve is reported in green (+50% α compared to the optimal operation line). Input parameter: T1 = 293 K, 0 ≤ α ≤ 73. The maximum value of α is the last one for which no numerical issues arise, due to the very high values of the temperatures and thus of the function f (T ). It is worth noting that realistic values of α range between 0 and 15.

Eq. (9) can be written by replacing T with θ, we define fˆ(T ) as fˆ(T ) =

Hθ , 1 − Kθ

fˆ0 (T ) =

H 2

(1 − Kθ)

(54)

Where H > 0, K > 0 are two constant parameters chosen to fit f (T ) within a range of temperature of interest, 0 < T < 343 K, resulting in H = 0.555 (dimensionless) and K = 0.011 1/K, with Tref = 273.15 K. Considering this function, the limit and critical values of µ are given by µ∗ =

(1 − Kθ3 )2 , H

µ∗ =

(1 − Kθ5 )2 , H

µc =

(1 − Kθ1 )(1 − Kθ6 ) H

Solving Eq. (52) and knowing that TH > T1 , TH takes the following expression when µ = µc : r α (1 − Kθ1 ) TH = T1 + K Since µc = 1/f 0 (TH ), this yields the following expression for µc : h i2 p −1 + Kθ1 + αK (1 − Kθ1 ) µc = H

(55)

(56)

(57)

which defines the optimal operation line of the system on the (µ, α) plane when varying α and T1 . The analytical solution intercepts both the µ and α axes at 2

µ0 =

(Kθ1 − 1) H 24

ACS Paragon Plus Environment

(58)

Page 25 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 14: Comparison between the actual and analytical solution determined by approximating f (T ) with fˆ(T ) with H = 0.555 and K = 0.011. (a) Comparison of functions f (T ) and fˆ(T ), 0 < T < 343 K. (b) Comparison of optimal operation curves for three values of T1 , 0 < α < 73; for the analytical solution the limit of α can be lower, to respect the constraint TH = T3 > T1 . Tref = 273.15 K.

and at α0 =

(Kθ1 − 1) K

(59)

The comparison between the actual and the analytical solutions is shown in Fig. 14. The two functions f (T ) and fˆ(T ) are compared in Fig. 14-a, whereas the resulting optimal operation curves for different T1 are compared in Fig. 14-b. One can note that the analytical form of µc describes well the actual solution, particularly within the range of α of interest. Indeed, the largest discrepancy is observed at very low values of α, where the system productivity approaches zero, and at very high values of α, which cannot be reached in practice and where the assumption of negligible distillate water flux is questionable. Furthermore, a better approximation is observed at higher values of T1 , which corresponds to smaller value of µ and larger values of productivity for a given value of α (see Section S4 of the Supporting Information). Overall, the analytical solution given by Eq. (57) provides an immediate tool for easily determining the optimal system operation for a given set of boundary conditions. 5. Conclusions This study investigates a desalination system based on a humidification-dehumidification (HDH) process driven by solar modules. The system is described within the framework of equilibrium theory by proposing a solution that is valid when assuming a constant saline water flow rate across the humidifier. This provides accurate predictions and valuable insight into physical processes, and allows understanding the optimal operation of the system, which is described by the air-to-water mass flow rate ratio µ, for given boundary conditions, which are specified by the solar energy source, α, and by the inlet temperature of the saline water, T1 . First, the equilibrium theory solution allows explaining the asymmetry between the humidifier and the dehumidifier, which does not result from mass transport limitations, but from the functional dependence of the enthalpy of moist air on temperature, which in turn depends on the functional dependence of the saturation pressure of air. On the one hand, humidifying air becomes easier when increasing the temperature, leading to simple wave interactions, to the corresponding smooth transitions, and to a greater entropy 25

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

generation (lower heat and mass exchange) in the humidifier. On the other hand, dehumidifying air becomes more difficult when decreasing the temperature, leading to shock interactions, to the corresponding sharp transitions, and to a smaller entropy generation in the dehumidifier. As a consequence, the minimum value of entropy generation is obtained in the dehumidifier when operating at a critical value of µ = µc ; here, both the outgoing air and water streams are in equilibrium with the corresponding incoming streams. Such a condition does not exist for the humidifier, which is characterized by two limit values of µ, namely µ∗ and µ∗ : below the former, the outgoing air is at equilibrium with the incoming water; above the latter, the outgoing water is at equilibrium with the incoming air. However, the simultaneous equilibrium at both ends of the exchanger cannot be achieved. This explains definitively the relevance of the dehumidifier in defining the optimal system operation. Then, the results of equilibrium theory are applied to the entire desalination process, where the dehumidifier and humidifier units are coupled through the solar modules. The optimal system operation is defined by determining the critical value µc for any set of boundary conditions α, T1 , i.e. by determining the value of µ that maximizes the production of clean desalinated water for given α and T1 . Such a maximum productivity curve separate a low-productivity and a high-productivity region, with a sudden drop in productivity occurring when going from the latter to the former. This behavior is well explained by equilibrium theory: for the coupled system, the critical value µc and the lower limit value µ∗ coincide, as they characterize the optimal operating points for dehumidifier and humidifier, respectively. This implies the sharp variation of the operation regime, hence of water productivity. Finally, equilibrium theory is used to determine an approximation of the optimal operating curve µc (α, T1 ), which can be calculated analytically. While no analytical solution is found when considering the actual exponential function describing the enthalpy of moist air, an analytical solution is found when using a rational function. The analytical form of µc describes well the actual solution, particularly within the range of α of interest, and provides an immediate tool for easily determining the optimal system operation. Notation ˜ ˆ

shock-related quantity approximate quantity

Symbols A cross sectional area [m2 ] solar area [m2 ] c specific heat capacity at constant pressure [kJ/(kgK)] f function of temperature characterizing moist air [K] G dry air mass flow rate [kg/s] H rational function parameter h latent heat of vaporization [kJ/kg] I solar irradiance [W/m2 ] K rational function parameter [1/K] L saline water mass flow rate [kg/s] p pressure [Pa] S˙ entropy generation [kJ/(kgK)] s˙ dimensionless entropy generation T temperature [K] t dimensional time coordinate [s] W evaporated/condensed water per unit dry air [kgw /kgda ] w moist air humidity ratio [kgw /kgda ] x dimensionless space coordinate z dimensional space coordinate [m]

26

ACS Paragon Plus Environment

Page 26 of 29

Page 27 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Greek letters α lumped parameter quantifying solar energy source [K]  liquid volume fraction η efficiency of solar modules λ propagation rate [m/s] µ air-to-water mass flow rate ratio ν air-to-water phase ratio π fresh water productivity ρ density [kg/m3 ] σ characteristic slope on physical plane [s/m] τ dimensionless time coordinate ω temperature derivative of function f Subscripts * lower limit ∞ state reached at infinite time + state just inside the exchanger, left-hand side state just inside the exchanger, right-hand side 0 initial A outlet air stream a inlet air stream B outlet liquid stream b inlet liquid stream c critical value D dehumidifier G dry air H humidifier L liquid ref reference U moist air v vapor Superscripts * per unit dry air upper limit Acronym CAOW closed-air open-water HCR heat capacity ratio HDH humidification-dehumidification PVT photovoltaic-thermal Supporting Information Supporting information is provided for further analysis concerning the comparison between the equilibrium theory solution and previous thermodynamic analyses (Section S1), the equilibrium theory solution of humidifier (Section S2), the relevance of the assumption of constant efficiency for the solar modules (Section S3), and the impact of inlet temperature of saline water. This information is available free of charge via the Internet at http://pubs.acs.org/.

27

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References [1] Bourouni, K.; Chaibi, M. T.; Tadrist, L. Water desalination by humidification and dehumidification of air: State of the art. Desalination 2001, 137, 167–176. [2] Narayan, G. P.; Sharqawy, M. H.; Summers, E. K.; Lienhard, J. H.; Zubair, S. M.; Antar, M. A. The potential of solardriven humidification–dehumidification desalination for small-scale decentralized water production. Renewable Sustainable Energy Rev. 2010, 14, 1187–1201. [3] Zamen, M.; Amidpour, M.; Soufari, S. M. Cost optimization of a solar humidification–dehumidification desalination unit using mathematical programming. Desalination 2009, 239, 92–99. [4] Abdelkareem, M. A.; El Haj Assad, M.; Sayed, E. T.; Soudan, B. Recent progress in the use of renewable energy sources to power water desalination plants. Desalination 2018, 435, 97–113. [5] M¨ uller-Holst, H.; Engelhardt, M.; Herve, M.; Sch¨ olkopf, W. Solarthermal seawater desalination systems for decentralised use. Renewable Energy 1998, 14, 311–318. [6] Elmutasim, S. M.; Ahmed, M.; Antar, M. A.; Gandhidasan, P.; Zubair, S. M. Design strategies of conventional and modified closed-air open-water humidification dehumidification systems. Desalination 2018, 435, 114–127. [7] Gabrielli, P.; Gazzani, M.; Novati, N.; Sutter, L.; Simonetti, R.; Molinaroli, L.; Manzolini, G.; Mazzotti, M. Combined water desalination and electricity generation through a humidification-dehumidification process integrated with photovoltaic-thermal modules: Design, performance analysis and techno-economic assessment. Energy Convers. Manage. X 2019, 1, 100004. [8] M¨ uller-Holst, H. Solar Desalination for the 21st Century: A Review of Modern Technologies and Researches on Desalination Coupled to Renewable Energies; Springer Netherlands: Dordrecht, 2007; pp 215–225. [9] Mistry, K. H.; Lienhard, J. H.; Zubair, S. M. Effect of entropy generation on the performance of humidification- dehumidification desalination cycles. Int. J. Therm. Sci. 2010, 49, 1837–1847. [10] Mistry, K. H.; Mitsos, A.; Lienhard V, J. H. Optimal operating conditions and configurations for humidification- dehumidification desalination cycles. Int. J. Therm. Sci. 2011, 50, 779–789. [11] Narayan, G. P.; Lienhard V, J. H.; Zubair, S. M. Entropy generation minimization of combined heat and mass transfer devices. Int. J. Therm. Sci. 2010, 49, 2057–2066. [12] Thiel, G. P.; Lienhard, J. H. Entropy generation in condensation in the presence of high concentrations of noncondensable gases. Int. J. Heat Mass Transfer 2012, 55, 5133–5147. [13] Chehayeb, K. M.; Narayan, G. P.; Zubair, S. M.; Lienhard, J. H. Thermodynamic balancing of a fixed-size two-stage humidification dehumidification desalination system. Desalination 2015, 369, 125–139. [14] Zamen, M.; Soufari, S. M.; Amidpour, M. Chemical Engineering Transactions; Italian Association of Chemical Engineering - AIDIC, 2011; Vol. 25. [15] Narayan, G. P.; Chehayeb, K. M.; McGovern, R. K.; Thiel, G. P.; Zubair, S. M.; Lienhard V, J. H. Thermodynamic balancing of the humidification dehumidification desalination system by mass extraction and injection. Int. J. Heat Mass Transfer 2013, 57, 756–770. [16] Narayan, P. G.; St. John, M. G.; Zubair, S. M.; Lienhard, J. H. Thermal design of the humidification dehumidification desalination system: An experimental investigation. Int. J. Heat Mass Transfer 2013, 58, 740–748. [17] Thiel, G. P.; Miller, J. A.; Zubair, S. M.; Lienhard, J. H. Effect of mass extractions and injections on the performance of a fixed-size humidification–dehumidification desalination system. Desalination 2013, 314, 50–58. [18] McGovern, R. K.; Thiel, G. P.; Prakash Narayan, G.; Zubair, S. M.; Lienhard V, J. H. Performance limits of zero and single extraction humidification-dehumidification desalination systems. Appl. Energy 2013, 102, 1081–1090. [19] Miller, J. A.; Lienhard V, J. H. Impact of extraction on a humidification-dehumidification desalination system. Desalination 2013, 313, 87–96. [20] Chehayeb, K. M.; Prakash Narayan, G.; Zubair, S. M.; Lienhard V, J. H. Use of multiple extractions and injections to thermodynamically balance the humidification dehumidification desalination system. Int. J. Heat Mass Transfer 2014, 68, 422–434. [21] Giwa, A.; Fath, H.; Hasan, S. W. Humidification-dehumidification desalination process driven by photovoltaic thermal energy recovery (PV-HDH) for small-scale sustainable water and power production. Desalination 2016, 377, 163–171. [22] Huang, X.; Li, Y.; Ke, T.; Ling, X.; Liu, W. Thermal investigation and performance analysis of a novel evaporation system based on a humidification-dehumidification process. Energy Convers. Manage. 2017, 147, 108–119. [23] Siddiqui, O. K.; Sharqawy, M. H.; Antar, M. A.; Zubair, S. M. Performance evaluation of variable pressure humidificationdehumidification systems. Desalination 2017, 409, 171–182. [24] Zubair, M. I.; Al-Sulaiman, F. A.; Antar, M.; Al-Dini, S. A.; Ibrahim, N. I. Performance and cost assessment of solar driven humidification dehumidification desalination system. Energy Convers. Manage. 2017, 132, 28–39. [25] Lawal, D. U.; Zubair, S. M.; Antar, M. A. Exergo-economic analysis of humidification-dehumidification (HDH) desalination systems driven by heat pump (HP). Desalination 2018, 443, 11–25. [26] Li, Y.; Huang, X.; Peng, H.; Ling, X.; Tu, S. Simulation and optimization of humidification-dehumidification evaporation system. Energy 2018, 145, 128–140. [27] Mazzotti, M.; Rajendran, A. Equilibrium Theory–Based Analysis of Nonlinear Waves in Separation Processes. Annu. Rev. Chem. Biomol. Eng. 2013, 4, 119–141. [28] Rhee, H.-K.; Aris, R.; Amundson, N. R. First-Order Partial Differential Equations, Volume I: Theory and Applications of Single Equations; Dover Publications, 1986; p 543. [29] Courant, R.; Friedrichs, K. O. Applied Mathematical Science 21 ; Springer: Berlin-Heidelberg-New York, 1979. [30] Whitham, G. B. Linear and Nonlinear Waves; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 1974.

28

ACS Paragon Plus Environment

Page 28 of 29

Page 29 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Intrinsic properties of humidificationdehumidification process

Optimal system operation Contour lines of water productivity

f’(T)>0 f’’(T)>0

Dehumidification Shock transition Humidification simple wave transition

29

ACS Paragon Plus Environment