Conjugate heat and mass transfer in a cross-flow hollow fiber membrane bundle used for seawater desalination considering air side turbulence

Author’s Accepted Manuscript Conjugate heat and mass transfer in a cross-flow hollow fiber membrane bundle used for seawater desalination considering air side turbulence Guo-Pei Li, Li-Zhi Zhang www.elsevier.com/locate/memsci

PII: DOI: Reference:

S0376-7388(17)30148-5 http://dx.doi.org/10.1016/j.memsci.2017.03.051 MEMSCI15162

To appear in: Journal of Membrane Science Received date: 17 January 2017 Revised date: 15 March 2017 Accepted date: 31 March 2017 Cite this article as: Guo-Pei Li and Li-Zhi Zhang, Conjugate heat and mass transfer in a cross-flow hollow fiber membrane bundle used for seawater desalination considering air side turbulence, Journal of Membrane Science, http://dx.doi.org/10.1016/j.memsci.2017.03.051 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Conjugate heat and mass transfer in a cross-flow hollow fiber membrane bundle used for seawater desalination considering air side turbulence Guo-Pei Li1, Li-Zhi Zhang1,2* 1

Key Laboratory of Enhanced Heat Transfer and Energy Conservation of Education Ministry, School of

Chemistry and Chemical Engineering, South China University of Technology, Guangzhou 510640, China. 2

State Key Laboratory of Subtropical Building Science, South China University of Technology, Guangzhou 510640, China.

*

Corresponding author. Tel/fax: 0086-20-87114268; Email: [email protected]

Abstract Fluid flow and heat and mass transfer in a hollow fiber membrane bundle for air humidification (MBH), which is the key process of a membrane desalination system, are investigated. The MBH is like a cross-flow shell-and-tube heat exchanger where heated seawater flows inside the fibers and the process air stream flows across the fiber bundle. Fluid flow and heat and mass transfer in the bundle are studied by considering the air side turbulence which exists even when the Reynolds numbers are as low as 50. Air side, water side and membrane side transfer equations are solved together in a conjugated way. Two turbulence models namely the standard k-ε model (STD k-ε) and the low-Reynolds-number k-ε model (Low Re k-ε) are adopted for the prediction of turbulent kinetic energy and its dissipation rate. In order to validate the numerical results and to compare

the

performance

of

the

different

turbulence

models,

a

membrane-based

air

humidification-dehumidification desalination (MHDD) system is designed and constructed to aid the investigation. Effects of module packing fractions and Reynolds numbers as well as the transverse and longitudinal pitches of the fibers on the mean Nusselt, Sherwood numbers and friction factors for air flow and water flow are analyzed. Correlations for performance analysis are regressed.

1

Graphical Abstract

H

Air in PL PD

A1

PT

A2

L do y

Solution in

z

W Computational cell

x

0.5 y*

Air side turbulent kinetic energy

0.25

0

0

0.5

1

1.5

2

x*

Keywords: Hollow fiber membrane humidifier; Desalination; Turbulence; Conjugate boundary condition; Heat and mass transfer

Nomenclature A

area (m2)

Am

fibers membrane area in the air stream side

Av

specific membrane area of unit volume of the module (m2/m3)

Atot

total membrane areas of the module (m2) 2

d

diameter (m)

dh

hydraulic diameter of the flow channels (m)

Dws

water diffusivity in saline water (m2/s)

Dva

moisture diffusivity in air (m2/s)

Dvm

moisture diffusivity in membrane (m2/s)

f

friction factor

H

module height (m)

Heva

evaporation heat of water (kJ/kg)

J

moisture emission rate (kg m-2 s-1)

L

fiber tube length (m)

k

turbulent kinetic energy (m2/s2)

nf

number of fibers

Nu

Nusselt number

p

pressure (Pa)

Pe

Peclet number

PD

diagonal pitch

PL

longitudinal pitch

PT

transverse pitch

Pr

Prandtl number

r

radius (m)

Re

Reynolds number

RH

Relative humidity (%)

Sc

Schmidt number

Sh

Sherwood number

T

temperature (K)

u

velocity (m/s)

W

module wide (m)

x/y/z X

coordinates for physical domain (m) mass fraction of water in saline water (kg water /kg saline water)

3

Greek letters δm

membrane thickness (m)

ε

dissipation rate of turbulent kinetic energy (m2/s3)

ζ*a

dimensionless temperature for periodic boundary condition

λ

heat conductivity (Wm-1K-1)

ν

kinematic viscosity (m2/s)

ξ*a

dimensionless humidity for periodic boundary condition

ρ

density (kg/m3)

φ

packing fraction

ψ

variable

ω

humidity (kg moisture/kg dry air)

Superscripts *

dimensionless

Subscripts a

air

b

bulk

c

cross section

e

equilibrium

eva

evaporation

h

heat

i

inlet

m

mass, membrane

max

maximum value

o

outlet

s

saline water

t

turbulence

tot

total

v

vapor 4

w

wall, water

x

x axis direction

y

y axis direction

1. Introduction Although our planet is often called the “Blue Planet”, warnings of increasing fresh water scarcity in the world are frequently heard. In fact, about 97.5% of water on the earth is sea or brackish water which contains large amounts of salt [1]. Meanwhile, the world is facing a global challenge to reliably supply its population with fresh water due to the shortages stemming from population growth, climate change, contamination of available freshwater supplies and public policy [2, 3]. As a solution, it suggests that the vast amounts of sea water could provide a nearly unlimited fresh water supply if coupled to an energy-efficient desalination technology. Of the various technologies for seawater desalination, the newly proposed membrane-based air humidification-dehumidification desalination (MHDD) is promising because its operating temperature can be less than 80°C, which would be powered by low grade energy sources such as solar energy and industrial waste heat, [4-6]. The technology relies on a humidification and a dehumidification process to get pure water from seawater. Compared to the dehumidification process, which has an efficiency as high as 0.94, the humidification efficiency is as low as 0.51 [7]. To intensify humidification efficiencies, hollow fiber membrane bundles for humidification (MBH, also called membrane contactors) are used to intensify moisture transfer. They have higher packing densities and thus higher heat and mass transfer capabilities [8-11]. According to their arrangement, the hollow fiber membrane bundles can be divided into parallel flow and cross flow modules. Cross-flow is commonly used because it’s convenient for duct sealing. Further, air side fluid flow and heat mass transfer are further intensified by the successive disturbances from the numerous fibers. As a result, the humidification effectiveness with this cross-flow arrangement is very encouraging [7, 11]. The bundle is like a cross-flow shell-and-tube heat exchanger. The heated saline water flows inside the fiber tubes, while the sweeping air flows across the fiber bank. The two fluids are separated 5

from each other by the semi-permeable membranes, which only selectively allow the permeation of water vapor but prohibit the transports of liquid solution through the membrane. The problems of liquid solution droplets cross-over and flow maldistribution in liquid side which are usually encountered in direct contact packed column humidification technologies, are thus overcome. Convective heat and mass transfer in hollow fiber membrane bundles are the key parameters for performance analysis of various membrane systems. Many researchers have investigated this subject experimentally and numerically. Either heat transfer [12, 13] or mass transfer [14, 15] in membrane contactors were studied for fuel cells or gas absorption applications. Others like Li et al. [16], Zhang et al. [17] and Huang et al. [18], considered the coupled heat and mass transfer through the porous membranes in the bundles for HVAC applications. However, as for the membrane bundles for desalination, information is still quite limited. Conjugate heat and mass transfer in membrane bundles for seawater desalination are different from the bundles in other applications by several features: (1) the operating air side Reynolds numbers are relatively high, i.e. Re>300, leading to turbulence inside the bundle. The turbulence is also contributed by the successive disturbances from the numerous fibers when the air flows across the fiber bank. Actually, the flow tends to be turbulent even in very lower Reynolds numbers [19, 20]. (2) The water side, the membrane, and the air side heat and mass transfer are closely coupled, so a uniform temperature or uniform heat flux boundary condition fails to exist. (3) Due to high packing fractions and inertial flow effects, the previously used free surface models are no longer valid, as it was indicated in [21]. Due to these reasons, the transport data built for other applications cannot be adopted for desalination uses. Turbulent fluid flow and heat and mass transfer in a membrane bundle for humidification in desalination applications will be the focus of this study. This is the novelty of this paper. Due to the less-development nature of turbulence under the Reynolds numbers around 300, the turbulence is transitional. Therefore two turbulence models namely the standard k-ε model (STD k-ε) and the low-Reynolds-number k-ε model (Low Re k-ε) are adopted to recognize and account for the degree of turbulence in bundle side. In addition to the continuity, momentum, energy and concentration equations, two additional equations, i.e., the turbulent kinetic energy and the dissipation rate equations are solved numerically to obtain the friction factors, Nusselt numbers and Sherwood numbers for air side. Those fundamental data is of interest for the effective design and optimization of the MHDD system. To validate the numerical results obtained and to compare the performance of the two different turbulence models, a 6

membrane-based air humidification-dehumidification desalination (MHDD) system is designed and constructed. Furthermore, the effects of transverse and longitudinal pitches of the fibers on the fluid flow and heat and mass transfer properties, as well as on the turbulent kinetic energy distribution in the bundle, are investigated.

2. Mathematical model 2.1. Governing equations The structure of the membrane bundle for humidification (MBH), where a numerous of hollow fiber tubes are packed orderly into a sealed shell, is like a shell-and-tube heat exchanger as shown in Fig. 1(a). Air stream flows in a cross flow manner across the hollow fiber tube bank. The physical structure and heat and mass transport data of the membrane fibers including membrane thickness (δm), fiber outer diameter (do), thermal conductivity of membrane (λm) and moisture diffusivity in membrane (Dvm) are listed in Table 1. Physical properties and inlet parameters for fluids are also listed. Due to the nature of periodicity in the air flow and the heat and mass transfer across the fiber bundle along the air flow direction, a typical periodic computational cell as shown in Fig. 1(b) is selected as the calculating domain. To simplify the calculation on the computational cell of complicate geometry with structured meshes, a boundary-fitted coordinate transforming method proposed in [22] is employed here to transform the irregular physical domain to a rectangular calculating domain as shown in Fig. 2(a) and Fig. 2 (b) respectively. Correspondingly, the mesh structures of the physical domains for both the solution side and air side are schematically plotted in Fig. 2 (c).

2.1.1 Solution side In the MBH, since the Reynolds numbers for the saline water stream in the fine fibers are generally less than 200 in practical applications, solution in the tube side can be considered to be laminar. Some other assumptions are made to simplify the mathematic model: it is fully developed for hydraulic, while temperature and concentration boundary layers are in developing along the channel length. In other words, the cross-sectional velocity field does not change, while temperature and concentration vary along with the flow direction. Based on the above assumptions, the governing equations comprising of the momentum equation, the energy equation and the mass equation are written as follows with dimensionless forms: 7

Momentum equation:

PL2  2us *  u2 s *  *2  *2  0 2 d h,s x y

(1)

Energy Equation: 2 * Ts d h,s us   2Ts  2Ts  =    zh* PL2us*  x*2 y*2 

(2)

2 * X s d h,s us   2 X s  2 X s  =  *2   zm* PL2us*  x*2 y 

(3)

x y and y*  PL PL

(4)

z z and zm*  Res Prs d h,s Res Scs d h,s

(5)

Mass Equation:

where

x* 

zh* 

where x*, y* and z* are dimensionless coordinates; Re, Pr and Sc are Reynolds, Prandtl and Schmidt numbers, respectively; dh,s is hydraulic diameter for the saline water flow channels (dh,s=di); PL and PT are the longitudinal and transverse pitches of the fiber bundle respectively. In Eqs. (1)-(3), us* and us* are the dimensionless local and mean velocities in the fiber tubes, T*s is the dimensionless temperature, where subscript “s” means saline water. X*s is the dimensionless mass fraction of saline water. They are defined as follows:

us*  

s sus d (dps dz )

us* 

Ts 

2 h,s

1 Ac

 u dA * s

Ts  Tai X  X ei and X s  s Tsi  Tai X si  X ei

where Ac is the cross-sectional area of the saline water flow channels (the fiber tubes), ν

(6)

(7)

(8)

is the kinematic

viscosity, p is pressure, Tai and Tsi are the temperatures of inlet air and inlet saline water respectively. Xsi is the inlet mass fraction of saline water (in kg pure water/kg saline water), and Xei is the equilibrium mass fraction

8

of the saline water with the air at the state of inlet temperature (Tai) and the inlet humidity (ωai). Those values used here have been listed in Table 1. NaCl solution with a mass fraction of 3.5% (kg NaCl/kg saline water) is used to substitute sea water in this research. The reason is that, on average, seawater in the world's oceans has a salinity of about 3.5% (kg salt/kg sea water). This means that every kilogram of seawater has approximately 35 grams of dissolved salts (predominantly NaCl). It is generally stable, so only the solution of 3.5% is investigated. The friction factor (fs):

fs  

d h,s dps dz

(9)

0.5s us 2

Reynolds number for saline water (Res):

Res 

d h,sus

(10)

s

For tube flow, the flow characteristic, convective heat transfer coefficients and convective mass transfer coefficients of the saline water can be represented by the product of the friction factor (fs) and the Reynolds number (Res), Nusselt numbers (Nu) and Sherwood numbers (Sh) respectively:

fs Res 

2

(11)

us 

Nus  

dTb,s 1 *   4 Tw,s  Tb,s  dzh,s

(12)



dX b,s 1 Shs   *   4  X w,s  X b,s  dzm,s

(13)

where Nus and Shs are the local Nusselt and Sherwood numbers for the saline water, subscripts “w” and “b” are “wall” and “bulk” respectively. Dimensionless bulk temperature and bulk mass fraction of saline water are obtained as follows:

Tb,s* 

* X b,s 

1  s

u Ac 1  s

u Ac

 u T dA

(14)

 u X

(15)

* * s s

* s

* s

dA

The mean values of the Nu and Sh over the entire fiber length are obtained simply as:

9

Nus

Shs

   

zh*

Nus dzh*

0



zh* 0

* zm

Shs dzm*

0



dz

* zm

0

dz

(16)

* h

(17)

* m

2.1.2 Air side The configuration of the problem under study is depicted in Fig. 1 (b). Though the air side Reynolds numbers are usually about 50-600 in practical applications, which is far below 2300, the critical transitional point from laminar to turbulent for round tube flow, it has been found that the flow tends to become turbulent due to the successive disturbances from the numerous fine fibers [19, 20]. A laminar model is not sufficient to disclose the detailed fluid flow and heat and mass characteristics in the bundle for air humidification. To address this problem, the standard k-ε model proposed by Jones and Launder [23], which has been proved to predict the more complex velocity fields better than other turbulence models for fiber bundles, is used. Due to its relative simplicity and adaptability but satisfactory performance, the STD k-ε model is still the most popular one used in numerical predictions of turbulence flow and heat transfer today [24, 25]. Although the model is a viable alternative to describe the role of the turbulent motions in high Reynolds-number flows, it is, nonetheless, difficult to formulate an exact predictions of turbulent convective heat and mass transfer immediately adjacent to a wall. Because there the local turbulent Reynolds numbers are very low, the direct influence of molecular viscosity on analyzing the velocity and the temperature fields in the near wall region cannot be ignored. In order to provide predictions of the flow within the viscous layer adjacent to the wall, Jones and Launder suggested that, the STD k-ε model must be enlarged. Subsequently, the low-Reynolds-number k-ε model (Low Re k-ε) was proposed in [26]. For comparison, both the STD k-ε model and the Low Re k-ε model are used. In addition, a laminar model is also used. The results of numerical calculation with the three models are plotted together for comparison. In the MBH, for the air side, flow and heat and mass transfer are uniformly distributed in the depth direction (z direction). Therefore, the flow is considered to be steady, incompressible, turbulent, and two-dimensional. The normalized governing equations of the fluid flow, heat and mass transfer for air stream are written as follows:

10

Continuity equation: * u*x u y  0 x* y*

(18)

Navier-Stokes equations: * * * (u*xu*x ) (uxu y )  pa*    a  t u*x     a  t u*x     a  t u y      2         x* y* 2x* x*  PLua x  y*  PLua y   y*  PLua x 

(19)

* *  pa*    a  t u y     a  t u y     a  t u*x    2       2y* x*  PLua x  y*  PLua y   x*  PLua y  

(20)

(u*xu*y ) x*



(u*y u*y ) y*



Energy and mass equations: * * (u*xTa* ) (u yTa )    a   t Ta*     a   t Ta*        x* y* x*  PLua x  y*  PLua y  

(u*xa* ) (u ya )   Dva  Dt a*    Dva  Dt a*        x* y* x*  PLua x  y*  PLua y   *

(21)

*

(22)

where, νa and αa are kinematic viscosity and thermal diffusivity of air stream respectively. Dva is the moisture diffusivity in the air. They can be found in Table 1. The turbulent dynamic viscosity νt and the turbulent thermal diffusivity αt as well as the turbulent moisture diffusivity Dt can be expressed as:

 t  cμ

k2



, t 

t Prt

and Dt 

t Sct

(23)

where Prt=0.95 and Sct =1.0. For STD k-ε model, the turbulent kinetic energy (k) and its dissipation rate (ε) appearing in Eq. (23) can be obtained by solving the k-equation and the ε-equation, which are put in the following way: *  (u*x k  ) (u y k )   k     k   ua t  PL i    Pk    k   k  x* y* x*  x  y*  y   PL ki ua ki *  (u*x  ) (u y ) ua t Pk  PL i  2                 c  c     ε ε 1 2 x* y* x*  x  y*  y   PL ki k  ua ki k 

(24)

(25)

where cμ=0.09, c1=1.44 and c2=1.92. However, for the Low Re k-ε model, the k-equation and the ε-equation are written as: *  (u*x k  ) (u y k )   k     k   ua t  PL i    Pk    D  k   k  x* y* x*  x  y*  y   PL ki ua ki

11

(26)

*  (u*x  ) (u y ) ua t Pk  PL i  2                 c  c  E  ε   ε  1 2 x* y* x*  x  y*  y   PL ki k  ua ki k 

(27)

where

 u* 2  u* y P  2  *x    *   x   y   k

2 a D  PLua 

u E  2 a t 3a PL  i 

  

2

  u* u*    *x  *y x   y 

  

2

(28)

  ( k  ) 1 2  2   ( k  ) 1 2  2       * *  x   y  

  2 u   2   2 u   2   2 u   *2x    *2x    *2y  x   y   x 

2

   2u y    *2   y

(29)

  

2

   

(30)

Compared with Eq. (24) for STD k-ε model, an extra term D* appearing in Eq. (26), which is equal to the dissipation rate in the immediate vicinity of the wall surface, is added to account for the fact that the dissipation processes are not isotropic [23, 26]. It is negligible in regions where the Reynolds number is high. Another extra term E* appearing in Eq. (27) is necessary to ensure the distribution of kinetic energy within the viscosity-affected region be in reasonable accordance with experiment [23, 26]. In addition, the terms cμ and c2 in Eqs. (23)-(27) which are constants in STD k-ε model, will become dependent upon the Reynolds number of turbulence (Ret) in the Low Re k-ε model:

cμ  0.09exp 2.5 / (1  Ret / 50)2 

(31)

c2  1.92 1  0.3exp( Ret2 ) 

(32)

where

Ret 

k2

(33)

 a

In Eqs. (24)-(27), the relevant diffusion coefficients, the dimensionless turbulent kinetic energy (k*) and the dimensionless turbulent dissipation rate (ε*) for the air stream are given [27]:

a  k 

k 

t

Prk and  ε  PLua

a 

t

Prε PLua

 ua d h,a  k 2k and      2 ki 0.01ua  i 500cμ ki2

where Prk=1.0 and Prε=1.3.

12

(34)

(35)

The hydraulic diameter of the air flow channels can be written as:

d h,a 

2( PT PL   ro2 )  ro

(36)

The packing fraction of the MBH:



 ro2

(37)

PT PL

In Eqs. (18)-(22), the dimensionless velocity, dimensionless temperature and dimensionless humidity of the air stream as well as the dimensionless pressure are defined as:

u *x 

Ta 

uy ux and u y*  ua ua

(38)

Ta  Tai   ai pa , a  a and pa*  Tsi  Tai si  ai 0.5a ua2

(39)

where ua is velocity of the air stream. It can be calculated from Rea= uadh,a/νa. ρa is air density, ωai is the humidity of inlet air; ωsi is the inlet air equilibrium humidity with saline water with a temperature Tsi and a mass fraction Xsi. They are listed in Table 1. Since the tube temperature and humidity are different with z direction, the conjugated air temperature and humidity on different x-y planes are different with z directions. The dimensionless bulk temperature and dimensionless bulk humidity at the inlet and outlet x-y planes are determined by:

 

* b,a

y

u *xTa*dy 

0

T



y

0

* a

u dy



and 

* b,a

 

y

u *xa*dy 

0



y

0

* a

u dy

(40)



The mean dimensionless bulk temperature and humidity for the whole flow are integrated along the entire z directions as follows:

* b,a

T

 

zh

0



* Tb,a dzh zh

0

dz

 h

* and b,a 



 zm

0



* b,a dzm  zm

0

dz

(41)

 m

Then the mean Nusselt and Sherwood numbers for the air stream are determined as [28]:   PT Rea Pra L Tb,ao  Tb,ai Nua   2Am Ta

13

(42)

Sha 

  PT Rea Sca L b,ao  b,ai  2Am a

(43)

where Am is the total area of membrane fibers in the air stream side, Am=2πroL. Subscripts “i” and “o” mean “inlet” and “outlet” respectively, T*a and ω*a are the log mean temperature and humidity differences between membrane surface and the air stream respectively, which can be expressed as:

Ta 

 Tao  Tai ao  ai  and    a  ln(Tao Tai ) ln(ao ai )

(44)

    and Tao  Tw,ao Tai  Tw,ai  Tb,ai  Tb,ao

(45)

     and ao ai  w,ai  b,ai  w,ao  b,ao

(46)

where subscripts “w” and “b” mean “wall” and “bulk”, respectively. T*ai and T*ao are the dimensionless temperature differences between the air stream and membrane surface at inlet and outlet planes. ω*ai and ω * ao

are the dimensionless humidity differences between the air stream and membrane surface at inlet and outlet

respectively. The mean friction factor (fa) used to predict the flow characteristic of air stream is given by [20]:

fa 

pa 2 N L (0.5a ua,max )

(47)

where NL is the number of fibers along the air stream flow direction; △pa is the air pressure drop between inlet and outlet. ua,max is maximum air stream flow rate in the tube bundle, it may occur at either the transverse plane (A1) or the diagonal plane (A2) [29], as indicated in Fig. 1(b):

PT ua , when PD   PT  do  2 PT  do

(48a)

PT ua , when PD   PT  do  2 2( PD  do )

(48b)

ua,max 

ua,max 

2.2. Boundary conditions The MBH studied has a staggered arrangement of cylindrical fiber tubes. The chosen periodic computational cell is shown in Fig. 1(b). Fig. 2(a) and Fig. 2 (b) show the irregular physical domain and the transformed rectangular calculating domain respectively. As shown, the computational cell is composed of three saline water flow channels (L1L2L3L4, M1M2M3M4 and N1N2N3N4) and an air flow channel 14

(ABCDEFGHI).Three membrane fiber tubes separate the two streams from each other. Planes IA and DE are respectively the entry and exit of air flow. Correspondingly, the mesh structures of the physical domains for both the solution side and air side are schematically plotted in Fig. 2 (c). The following are the boundary conditions necessary to complete the numerical computation. No-slip boundary conditions are imposed on both the inner surface and outer surface of membrane: At planes AB, CD, FGH, L3L4, M3M4 and N3N4: u*x  0 , uy*  0 and us*  0

(49)

Symmetric conditions of flow, temperature, humidity and saline water mass fractions are imposed at planes L4L1, L1L2L3, M4M1M2, M2M3, N4N1N2N3 for the saline water and planes BC, EF, HI for the air stream. Therefore the gradient of the flow parameters at the normal direction are equal to zero:

 * * * *  0 , ψ= u , T , ω a and X s n

(50)

where n refers to the normal direction. Considering the nature of periodicity of the air flow channel, periodic inflow and outflow boundaries of velocity are:

u 

* x IA

  u*x 

DE

and

u 

* y IA

  uy* 

DE

(51)

In the fiber bundles, after flowing through 4-6 rows of fiber tubes, the air stream is fully developed both in velocity and thermally [30] while there are 34 rows of fiber tubes in the MBH in this research. Actually in practical application, for the whole hollow fiber membrane bundle for air humidification, we are more concerned about the average heat and mass transfer coefficients. So the assumption that the air flow is fully developed before it enters the periodic computational cell is made to simplify the mathematic model. The definition of the thermally and concentration fully developed regime is that the shapes of the dimensional temperature and humidity profiles rather than the specific values at successive streamwise locations are the same. When the heat and mass transfer boundary layers are fully developed, the dimensionless temperature and humidity defined below are employed for the periodic boundary condition. The defined values on planes IA and DE are equal as:

 

* a IA

  a* 

DE

and

 

* a IA

 a* 

DE

(52)

where ζ*a and ξ*a are dimensionless temperature and humidity of the air employed for the periodic boundary condition, which are given by [28, 31]:

15

 a 

* Ta*  Tw,a * * Tb,a  Tw,a

and a 

* a*  w,a * * b,a  w,a

(53)

where subscripts “w” and “b” mean “wall” and “bulk”, respectively. T *w,a and ω *w,a are dimensionless temperature and humidity of the air at the wall (membrane) surface respectively; T *b,a and ω *b,a are dimensionless bulk temperature and dimensionless bulk humidity of the air stream respectively. The inlet boundary conditions for saline water are:

Ts*

* zh,s 0

 1 and X s*

* zm,s 0

1

(54)

The inlet boundary conditions for air are [27]:

Ta* k

 x  0

x  0

 0 and a*

0.01ur2 2 and    ki

x  0

x  0

0



(55)

500cμ k 2 ua d h,a  i

(56)

where ur is resultant velocity, the ki and εi can be calculated by Eq. (35). It should be clarified that the boundary conditions of the humidification process in the hollow fiber membrane bundle for seawater desalination are neither uniform heat flux nor uniform temperature. It is naturally-formed boundary conditions under which the heat and mass transfer are strongly coupled because of water vapor evaporation on membrane surface inside the fiber tubes. In other words, the saline water side heat flux is the sum of the air side sensible heat flux and the latent heat flux. Since the membranes are rather thin and the heat conductivity is relatively high, trans-membrane temperature differences are negligible [32]. Therefore the normalized heat balance equation and the mass balance equation on the membrane surfaces are:

Ts* n

= m,s

a Ta*  s n

*  H eva  m,a

a* n

J m,s =J m,a  J m

(57) m,a

(58)

where n refers to the normal direction, λa and λs are the heat conductivity of the air and the saline water respectively. H*eva is the dimensionless evaporation heat. Jm,s and Jm,a are the moisture emission of the saline water side and air side on membrane surface respectively, Jm is the moisture emission through the membrane. They are defined as： * H eva =

a Dva H eva  si  ai    s  Tsi  Tai  16

(59)

J m,s   s Dws

X s n

and J m,a   a Dva m,s

J m  a Dvm

m,s  m,a m

a n

(60) m,a

(61)

where Dva is moisture diffusivity in air, Heva is the evaporation heat of water (2501 kJ/kg). Dws and Dvm are the water diffusivity in saline water and the moisture diffusivity in membrane respectively, δm is the membrane thickness, ωm,s and ωm,a are the humidity of the saline water side and air side on membrane surface respectively.

2.3. Numerical solution procedure The normalized governing partial differential equations, together with the boundary conditions of the mathematical model for turbulent flow, are given in sections 2.1 and 2.2 respectively. In this section, the detailed solution procedure will be introduced. In order to obtain the numerical solution of the mathematical model, as a first step, a body-fitted coordinate system is establish to transform the irregular physical domain to the rectangular calculating domain. Then those governing equations as well as the calculating domain are discretized in the structured meshes by means of a finite volume method. Considering that heat and mass transfer are closely coupled through the membrane, and the velocity field is related to the temperature and concentration fields, an alternate direction implicit (ADI) method is used to solve these discretization equations. The equations for the air stream are solved at each x-y plane across the z direction. Fig. 3 demonstrates the detailed flow chart of the numerical calculation. A grid independence test is done to ascertain that the discretization error is reduced to an acceptable level. It can be found that 21×81 grids for air flow channel (ABCDEFGHI), 21×21 grids for the left and the right saline water flow channel (L1L2L3L4 and M1M2M3M4), 21×31 for the middle saline water flow channel (N1N2N3N4) and 41 grids in the axial direction had enough numerical accuracy (less than 0.1%) compared with 31×101 grids for air flow, 31×31 grids for the left and the right saline water flow, 31×41 for the middle saline water flow and 61 grids in the axial direction.

3. Experimental study For the sake of verifying the validity of the present computational model, a membrane-based air 17

humidification-dehumidification desalination (MHDD) experimental system is designed and constructed. As shown in Fig.4, the whole test rig consists of a heat storage water tank with an electric heater, a MBH and a dehumidifier (a fin-and-tube heat exchanger). The outside of the heat storage tank is covered by a 5 cm thick polyurethane foam to insulate heat transfer from the inside to the environment. The inner shell of the tank is made of 304 stainless steel to ensure the equipment is corrosion resistant to saline water. The MBH proposed in this research is like a shell-and-tube heat mass exchanger, where solution flows in tube side and air flows across the fiber bundle. NaCl solution with a mass fraction of 3.5% is used to substitute sea water in the experiment. The saline water stored in the heat storage tank is heated to 70 °C (Tsi) by the electric heater and it then flows through the fibers, while the fresh air at 30 °C (Tai) powered by a centrifugal fan with frequency conversion governor flows outside the fibers in the shell side of the module. The fibers are made with porous polymer PVDF (polyvinylidene fluoride). The hydrophobic PVDF layer mainly provides the mechanical strength. It also prevents liquid molecules from entering the pores of membranes which has been verified by a long time experiment. The system is based on an open cycle for air stream and a closed cycle for saline water. Under water vapor pressure differences between the two sides of the membrane, moisture is driven from the water through the membrane pores and eventually reaches the fresh air stream. The heated and humidified air from the humidifier is then cooled and dehumidified by the dehumidifier; meanwhile liquid fresh water is produced. The saline water leaving the MBH is collected in the water tank and it is then recalculated to the MBH after heated again. To ensure the saline water has a constant concentration during the test, the produced fresh water collected in the fresh water storage tank flows back into the heat storage water tank. The MBH and the saline water lines are covered by PVC (polyvinyl chloride) heat insulation cotton to reduce the heat losses. Several screens are placed along the air duct to uniformly distribute the air stream. A Data Transmission System (DTS) is used to collect temperature and humidity data of the fluids and then transfer them to a computer. In addition, in order to maintain constant inlet air states, the whole system is installed in an air-conditioned chamber. To compare with the theoretical results, several key parameters of the system are needed to be measured, such as the temperature of saline water at inlet and outlet of the humidifier, the relative humidity and temperature of air at inlet and outlet of the humidifier. Therefore two PT100 temperature sensors with an accuracy of ±0.1°C are used to measure the saline water temperature at inlet and outlet of the humidifier. To measure the flow rates of water entering the humidifier, a rotor-flow meter is used with an accuracy of ±0.5%. 18

Three digital temperature-humidity sensors (AF3485Y) with an accuracy of ±2.5% and a resolution of 0.1% are employed to measure the temperature and relative humidity of the air. Air flow rates are obtained by a hot wire anemometer (TESTO 425) with an accuracy of ±0.03m/s. Tube side and shell side pressure drops are measured by electronic pressure gauges with an accuracy of ±0.1Pa. All sensors are calibrated before the experiment. A conductivity meter (DDS-307) with accuracy of ±0.5% combined with a DJS-10C platinum black electrode is used to measure the electrical conductivity of the saline water at inlet and outlet of the humidifier. The mass fractions of the solution to and out the membrane bundle are obtained by titration method using silver nitrate solution. The conductivity and the concentration are all measured directly in our laboratory. Based on those measured values, then a linear equation is regressed for the NaCl solution at 28 °C, as described in [7]. Then the mass fraction of the saline water can be obtained from the measured electrical conductivity using the experimental correlation conveniently. With those measured inlet and outlet parameters of the two streams, the mean Nusselt and Sherwood numbers and the friction factor can be estimated by Eqs. (42)-(44) and (47), with temperature and humidity difference substituted by the log mean temperature and humidity difference between the inlets and the outlets respectively. Based on the accuracy of each measuring apparatus, the final uncertainties of the experimental results are estimated according to Ref. [33]:

Δx Δy f f 2 Δx2 2 f 2 Δxn 2 2  [( ) 2 ( 1 ) 2  ( ) ( )  ( ) ( ) ] y x1 y x2 y xn y 1

(62)

where f is the function of the independent variables, x1, x2, etc. is the variables, x is the absolute error associated with the variables, y y stands for the relative error. Therefore the uncertainty for the concentration of the salt solution estimated by the measured electrical conductivity is 0.51%, and the final uncertainties for the tested friction factor, Nusselt and Sherwood numbers are 3.6%, 7.5% and 7.9%, respectively.

4. Results and discussion 4.1. Model validations In order to verify this model, a comparison study between the theoretical results and experimental data under the naturally-formed boundary conditions is performed for both heat mass transfer and flow resistance 19

in tube banks. The MBH used in the test rig has a staggered arrangement of cylindrical fiber tubes with equal transverse pitch (PT) and longitudinal pitch (PL) of PL/do= PT/do=1.75, and the module packing fractions (φ) is 0.256. Other designed operating conditions are listed in Table 1. The saline water in the heat storage water tank is heated firstly to the required temperature by the electric heater in the tank and then the hot saline water is transported to the tube side of the humidifier by a corrosion resistant magnetic circulating pump. During the test, the air outlet temperature (Tao), relative humidity (RHao) and pressure drop (△Pa) in the MBH under designed operating conditions are measured. Based on those measured value, then the experimental values of the mean Nusselt numbers ( Nua ) and Sherwood numbers ( Sha ) and the friction factor (fa) can be estimated based on Eqs. (42)-(44) and (47). Comparisons are made between the calculated data from STD k-ε model and from Low Re k-ε model as well as the experimentally results under different air Reynolds numbers ranging from 50 to 600 which is common in practical operations. For comparison, besides the proposed two turbulence models, a laminar model is also used. They are plotted in Figs. 5, 6 and 7 respectively. In general, either for the Nusselt numbers, the Sherwood numbers or for the friction factors, the deviations of the calculated results with different mathematical models increase gradually with increasing Rea. As seen, the Low Re k-ε model shows a better agreement with the experiments in the tested Rea ranges from 50 to 600, especially when Rea is greater than 300. The laminar model only fits the tests well at Reynolds numbers below 300. However, a noticeable deviation is observed between the results of STD k-ε model and the experiments after Rea is higher than 100. For the Low Re k-ε model, the deviations between the experimental results and the numerical data are within 4.9% for Nusselt numbers, 4.6% for Sherwood numbers and 5.9% for friction factor. In comparison, the deviations of the three parameters with the STD k-ε model are 18.5%, 15.9% and 14.3% respectively. In Figs.8 and 9 are plotted respectively the contours of the dimensionless turbulent kinetic energy (k*) and its dissipation rate (ε*) in the periodic computational cell with a staggered arrangement of PT/do=PL/do=1.5 and an air Reynolds number of 400, by comparing the predictions with the two turbulence models. The velocity vectors for the air are also plotted in Fig. 8 (a). It is seen that the distributions of k* and ε* are similar to the velocity fields. The values of k* and ε* are low in the upstream and downstream region of the tubes where the velocities are low. Large values of k* and ε* are found in the channel center and in the flow separation regions where velocities are higher. Moreover, compared to the predictions by the Low Re k-ε model, the values of k* and ε* predicted by the STD k-ε model are much greater. It means that the turbulence intensities predicted by

20

the STD k-ε model are higher than those by the Low Re k-ε model, which may have contributed to the 20% overestimation of the Nusselt numbers, Sherwood numbers and friction factor. In addition, the calculated performance data are compared to those from references, as given in Table 2. The pitch-to-diameter ratio is PT/do=PL/do=1.5. The calculated friction factors and Sherwood numbers that are under naturally formed boundary conditions can be compared directly with reference [28, 34, 35]. References [28, 34] provide only the Nusselt numbers under uniform temperature and uniform heat flux boundary conditions, so the Nusselt numbers under uniform temperature and uniform heat flux boundary conditions are deliberately calculated by the model developed here, but the boundary conditions are changed from naturally formed (conjugated) to uniform temperature or heat flux. It can be found that, the mean Nusselt numbers under uniform temperature and uniform heat flux boundary conditions in this study are large than those from references under the same thermal boundary conditions. The friction factors predicted by the proposed Low Re k-ε model are also greater than [28, 34]. The maximum difference is within 4.8% for ( Nua )T, 1.7% for ( Nua )H and 6.4% for fa respectively. The deviations may contributed by the turbulence in the bundle which has intensified fluid flow and heat and mass transfer. The comparison indicates that the mathematical model is successful in predicting the flow and heat transfer in the cross flow MBH. It can also be seen that the current Sherwood numbers are higher than those for pure air dehumidification applications that use 35% LiCl solution as liquid desiccant [35]. The reason is that, the concentration entry length for the solution flow in the circular inner tubes is 29.8cm in this application, which is longer than that for LiCl solution (25cm) in [35]. For heat and mass transfer under naturally formed boundary conditions, the thermal and concentration entry lengths of the solution stream in the tube side have a great effect on the air side Nusselt numbers and Sherwood numbers. The detailed analysis of the tube side can be seen in section 4.5.

4.2. Air side: effects of the module packing fractions and Reynolds numbers The flow and heat mass transfer characteristics in the tube bundle, which are of interest in the practical applications, are investigated in this section. Considering that the Low Re k-ε model has a satisfactory agreement with the present experimental results compared to the laminar model and the STD k-ε model, when covering a wider range of Reynolds numbers, the following sections use this model to investigate the performances of the MBH under various module packing fractions and operating conditions. The variations of mean Nusselt numbers and mean Sherwood numbers for the air stream with Rea ranging 21

from 50-600 at different module packing fractions are illustrated in Figs. 10 and 11 respectively. The range of module packing fraction is from 0.126 to 0.503, which is typically employed in practical applications. It is clear that the Nua always increases with Rea with different packing fractions. Similar trends are observed for the air side Sha . This is due to the fact that the higher the values of Re is, the greater the turbulence intensities are. Higher Res make the flow more complicated and mixed. Heat and mass transfer in the tube bundle occur primarily by bulk flow convections. Therefore the convective heat and mass transfer in the bundle are enhanced and thus the mean Nusselt numbers and Sherwood numbers become larger. However this increase has different slopes at different Res. The effect of Reynolds numbers on the mean Nusselt numbers and mean Sherwood numbers are more significant for more compact bundles. The effects of Reynolds numbers and module packing fractions on the friction factors (fa) are indicated in Fig.12. As seen, the friction factor decreases steeply when Rea ranges from 50 to 250. Then the decrease becomes slower after the Rea is higher than 300. Moreover, the packing fraction of the modules also has a great effect on the friction factor. Increasing the packing fraction leads to an increase in the friction factor. This behavior may be attributed to the sharp reduction in the hydraulic diameter of the air flow channel, which would increase the turbulence intensities in the narrowed spaces.

4.3. Air side: effects of transverse and longitudinal pitches The effects of the transverse pitch (PT) and longitudinal pitch (PL) of the fibers on heat and mass transfer and friction factors are shown in Figs. 13, 14 and 15. As seen in Fig. 1(b), PT is the transverse tube pitch perpendicular to the direction of flow. PL is the longitudinal pitch along the air flow direction. As shown in Figs. 13, 14 and 15, heat and mass transfer increase mainly with decreasing transverse pitches, but to a lesser extent with decreasing longitudinal pitches. Similarly, compared to the longitudinal pitch, the transverse pitch has a more significant impact on friction factors. To disclose the different working-mechanisms of transverse and longitudinal pitch, the typical profiles of dimensionless turbulent kinetic energy (k*) in the air streams for different transverse pitches (PT) and longitudinal pitches (PL) are plotted for comparison, as shown in Fig. 16. As seen, decreasing the transverse pitch is more influential than decreasing the longitudinal pitch in having an increased turbulent kinetic energy. It means that the turbulence intensity is more intense in the tube bundle with small transverse pitches, which consequently results in a better enhancement in the mean Nusselt and Sherwood numbers. 22

4.4. Air side: correlations for performance evaluation In this section, in order to make it convenient to predict the heat and mass transfer as well as the pressure drop of the bundle in the design and optimization of the MBH, which is a key component in the membrane-based air humidification-dehumidification desalination system, a set of correlations are regressed based on the above mentioned numerical results under naturally formed boundary conditions. The simple calculation of mean Nusselt numbers, Sherwood numbers and friction factor can be written as:

Nua  CRem Pr1 3

(63)

Sha  DRen Sc1 3

(64)

fa  a  Reb

(65)

In Eqs. (63)-(65) the parameters C, D, m, n, a and b are constants depended on the geometrical configurations. They are related to the packing fractions of the module which are given in tabular forms in Table 3. It is needed to specify that those correlations are only available for air stream in tube bundle with staggered arrangement of fibers (PT= PL). The application range is for the module packing fractions ranging from 0.126 to 0.503 and for Reynolds numbers ranging from 50 to 600. The data for inline arrangements of fibers (PT=PL) are given in Table 4. Compared to staggered arrangements, MBH with inline arrangements usually have lower friction factors and lower heat and mass transfer coefficients .The maximum deviations of these correlations with the calculated values are within 3.6% for Nua , 3.1% for Sha and 9.3% for fa respectively.

4.5. Tube side laminar flow As mentioned previously, the flow and heat mass transfer characteristics in the tube bundle are of interest in the practical applications and are investigated in detail above. In this section, the tube side transport characteristics are also analyzed. For the saline water stream in the fine fibers, since the Reynolds numbers are generally less than 200 in practical applications, solution in the tube side can be considered to be laminar. The calculated value of fsRes (from Eq. (11)) is 64.68, which is consistent with the previous numerical values of 64 from [34]. The fully developed local Nusselt numbers of the saline water stream in circular ducts (from Eq. (12)) are 3.74 and 4.43 under uniform temperature and uniform heat flux boundary conditions respectively,

23

which are also consistent with the values of 3.66 and 4.36 from [34]. The local Nusselt numbers under naturally formed boundary conditions decrease sharply from large values at the inlet to fully developed values after the entry length. The fully developed mean Nusselt numbers and Sherwood numbers are 4.71 and 5.48 in this study which are larger than the values of 4.56 and 5.32 in air dehumidification application which use 35% LiCl solution as liquid desiccant in [35]. The reason is that the length of the thermal inlet region and the concentration inlet region are z=2.63cm (z*h =0.245) and z=29.8cm (z*m=0.051) for the saline water, which are larger than the values of 0.5cm and 25cm in [35]. As seen, in contrast to the rapid development of the thermal boundary layer, the concentration boundary layer develops much slower. The reason is that the Schmidt number of NaCl solution is relatively larger (Scs=150) while the Prandtl number is only 2.75. However the length of the entire fiber tube is 35cm which is long enough for the saline water stream to get thermal and concentration fully developed.

5. Conclusions A hollow fiber membrane bundle for humidification (MBH) is used for seawater desalination. The air side convective heat and mass transfer as well as the fluid flow characteristics are investigated considering turbulence in the bundle. Two turbulence models namely the standard k-ε model (STD k-ε) and the low-Reynolds-number k-ε model (Low Re k-ε) are adopted for the prediction of the turbulent motions. A membrane-based air humidification-dehumidification desalination (MHDD) system is also designed and constructed to validate the proposed mathematical model. Following conclusions can be given: (1) Generally, the Low Re k-ε model shows a better agreement with the experiments than the laminar and the STD k-ε model in the tested Rea range from 50 to 600, especially when Rea is greater than 300. The laminar model only fits the tests well at Reynolds numbers below 300. A noticeable deviation is observed between the STD k-ε model and experiments after the Rea is higher than 100. For the Low Re k-ε model, the deviations between the experimental results and the numerical data are within 6.0%, while the deviations are about 20% for the STD k-ε model. This is due to the reason that the in-bank turbulent kinetic energy and its dissipation rate are exaggerated by the STD k-ε model. (2) Tube side: in contrast to the rapid development of the thermal boundary layer, the concentration boundary layer develops much slower. (3) Air side: Reynolds numbers and the module packing fractions have a great impact on the fluid flow and heat mass transfer performance. With increases in Reynolds numbers and packing fractions, the 24

convective heat and mass transfer in the bundle are enhanced, while the frictions are also increased. (4) Air side: Heat and mass transfer increase mainly by decreasing transverse pitches, rather than by decreasing longitudinal pitches. Shortening the transverse pitch has a more significant impact on the turbulent kinetic energy distribution and the turbulence intensity compared to shortening the longitudinal pitches. (5) Turbulence should be considered when modeling the membrane bundle.

Acknowledgement This Project was supported by the National Science Fund for Distinguished Young Scholars of China, No. 51425601. It was also supported by Natural Science Foundation of China, No. 51376064.

References [1] T. Oki, S. Kanae, Global hydrological cycles and world water resources, Science 313 (2006) 1068-1072. [2] C.J. Vörösmarty, P. Green, J. Salisbury, R.B. Lammers, Global water resources: vulnerability from climate change and population growth, Science 289 (2000) 284-288. [3] K.N. Knust, D. Hlushkou, U. Tallarek, R.M. Crooks, Electrochemical Desalination for a Sustainable Water Future, ChemElectroChem 1 (2014) 850-857. [4] T. C. Chen, C. D. Ho, Immediate assisted solar direct contact membrane distillation in saline water desalination, J. Membr. Sci. 358(2010) 122-130. [5] E. Guillén, J. Blanco, G. Zaragoza, et al. Experimental analysis of an air gap membrane distillation solar desalination pilot system, J. Membr. Sci. 379(2011) 386-396. [6] L.Z. Zhang, G.P. Li, Energy and economic analysis of a hollow fiber membrane-based desalination system driven by solar energy, Desalination 404 (2017) 200-214. [7] G.P. Li, L.Z. Zhang, Investigation of a solar energy driven and hollow fiber membrane-based humidification–dehumidification desalination system, Appl. Energy 177 (2016) 393-408. 25

[8] L.Z. Zhang, S.M. Huang, Coupled heat and mass transfer in a counter flow hollow fiber membrane module for air humidification, Int. J. Heat Mass Transfer 54 (2011) 1055-1063. [9] S. Al-Obaidani, E. Curcio, F. Macedonio, G. Di Profio, H. Al-Hinai, E. Drioli, Potential of membrane distillation in seawater desalination: thermal efficiency, sensitivity study and cost estimation, J. Membr. Sci. 323 (2008) 85-98. [10] K. Kneifel, S. Nowak, W. Albrecht, R. Hilke, R. Just, K.V. Peinemann, Hollow fiber membrane contactor for air humidity control: Modules and membranes, J. Membr. Sci. 276 (2006) 241-251. [11] L.Z. Zhang, Coupled heat and mass transfer in an application-scale cross-flow hollow fiber membrane module for air humidification, Int. J. Heat Mass Transfer 55 (2012) 5861-5869. [12] N. Djilali, D. Lu, Influence of heat transfer on gas and water transport in fuel cells, Int. J. Therm. Sci. 41 (2002) 29–40. [13] T. V. Nguyen, R.E. White, A water and heat management model for Proton-Exchange-Membrane fuel cells, J. Electrochem. Soc. 140 (1993) 2178-2186. [14] S.R. Wickramasinghe, M.J. Semmens, E.L. Cussler, Mass transfer in various hollow fiber geometries, J. Membr. Sci. 69 (1992) 235-250. [15] S. Shirazian, A. Moghadassi, S. Moradi, Numerical simulation of mass transfer in gas–liquid hollow fiber membrane contactors for laminar flow conditions, Simul. Model. Pract. Th. 17 (2009) 708-718. [16] T. Li, N.G. Deen, J.A.M. Kuipers, Numerical investigation of hydrodynamics and mass transfer for in-line fiber arrays in laminar cross-flow at low Reynolds numbers, Chem. Eng. Sci. 60 (2005) 1837-1847. [17] L.Z. Zhang, S.M. Huang, L.X. Pei, Conjugate heat and mass transfer in a cross-flow hollow fiber membrane contactor for liquid desiccant air dehumidification, Int. J. Heat Mass Transfer 55 (2012) 8061-8072. [18] S.M. Huang, L.Z. Zhang, L.X. Pei, Transport Phenomena in a Cross-Flow Hollow Fibre Membrane Bundle Used for Liquid Desiccant Air Dehumidification, Indoor Built Environ. 22 (2013) 559-574. [19] L. Zhang, W. Zhong, J. Chen, J. Zhou, Fluid Flow and Heat Transfer in Plate-Fin and Tube Heat Exchangers in a Transitional Flow Regime, Numer. Heat Transfer Appl. 60 (2011) 766-784. [20] L.Z. Zhang, Z.X. Li, Convective mass transfer and pressure drop correlations for cross-flow structured hollow fiber membrane bundles under low Reynolds numbers but with turbulent flow behaviors, J. Membr. Sci. 434 (2013) 65-73. 26

[21] V.A. Kirsch, V.I. Roldugin, A.V. Bildukevich, V.V. Volkov, Simulation of convective-diffusional processes in hollow fiber membrane contactors, Separation & Purification Technology, 167 (2016) 63-69. [22] J.L. Niu, L.Z. Zhang, Heat transfer and friction coefficients in corrugated ducts confined by sinusoidal and arc curves, Int. J. Heat Mass Transfer 45 (2002) 571-578. [23] W.P. Jones, B.E. Launder, The prediction of laminarization with a two-equation model of turbulence, Int. J. Heat Mass Transfer 15 (1972) 301-314. [24] W.E. A Balabel, On the performance of linear and nonlinear turbulence models in various jet flow applications, Eur. J. Mech. B/Fluids 30 (2011) 325-340. [25] F. Juretić, H. Kozmar, Computational modeling of the neutrally stratified atmospheric boundary layer flow using the standard k–ε turbulence model, J. Wind Eng. Ind. Aerod. 115 (2013) 112-120. [26] W.P. Jones, B.E. Launder, The calculation of low-Reynolds-number phenomena with a two-equation model of turbulence, Int. J. Heat Mass Transfer 16 (1973) 1119-1130. [27] S.V. Patankar, E.M. Sparrow, M. Ivanović, Thermal interactions among the confining walls of a turbulent recirculating flow, Int. J. Heat Mass Transfer 21 (1978) 269-274. [28] W.M. Kays, M.E. Crawford, Convective heat and mass transfer, McGraw-Hill, Inc, New York, 1993. [29] T.L. Bergman, F.P. Incropera, A.S. Lavine, Fundamentals of heat and mass transfer, John Wiley & Sons, 2011. [30] A. Žukauskas, Heat transfer from tubes in crossflow, Advances in heat transfer, 8 (1972) 93-160. [31] S.V. Patankar, C.H. Liu, E.M. Sparrow, Fully developed flow and heat transfer in ducts having streamwise-periodic variations of cross-sectional area, J. Heat Transfer, 99 (1977) 180-186. [32] L.Z. Zhang, Heat and mass transfer in a cross-flow membrane-based enthalpy exchanger under naturally formed boundary conditions, Int. J. Heat Mass Transfer 50 (2007) 151-162. [33] I. Guttman, S.S. Wilks, H.J. Stuart, Introductory Engineering Statistics, John Wiley & Sons, New York, 1965. [34] F.P. Incropera, D.P. Dewitt, Introduction to Heat Transfer, third ed., John Wiley & Sons Inc, New York, 1996. [35] S.M. Huang, L.Z. Zhang, L.X. Pei, Transport Phenomena in a Cross-Flow Hollow Fibre Membrane Bundle Used for Liquid Desiccant Air Dehumidification, Indoor Built Environ. 22 (2013) 559-574.

27

Table 1 Physical and transport properties of the membrane bundle under design operating conditions. PL/do= PT/do=1.75, φ=0.256.

Symbol

Unit

Value

Symbol

Unit

Value

W

mm

90

Dvm

m2/s

9.0×10−7

H

mm

190

Tai

°C

30.0

L

mm

350

Tsi

°C

70.0

di

mm

1.3

Res

-

30

do

mm

1.5

ωai

kg/kg

0.0174

δm

mm

0.1

ωsi

kg/kg

0.2692

2477

Xsi

kg/kg

0.965

nf

2

3

Av

m /m

682.6

Xei

kg/kg

0.74

Atot

m2

4.08

νa

m2/s

1.589×10−5

ρa

kg/m3

1.16

νs

m2/s

4.5×10−7

ρs

kg/m3

977.5

Pra

-

0.71

λa

Wm−1K−1

0.0263

Prs

-

2.75

λs

Wm−1K−1

0.64

Sca

-

0.564

Dws

2

m /s

Scs

-

150

Dva

2

3.0×10

m /s

−9

2.82×10

−5

Table 2 Comparisons of mean Nusselt numbers under uniform temperature condition ( Nua )T and under uniform heat flux boundary conditions ( Nua )H, and the friction factors (fa) as well as the mean Sherwood numbers ( Sha ) under naturally formed boundary conditions for the air flow across the tube banks, PT= PL, Rea=100. ( Nua )T φ

Sha

fa

( Nua )H

PL/do (This

Refs.

(This

Refs.

(This

Refs.

(This

Ref.

study)

[28, 31]

study)

[28, 31]

study)

[28, 31]

study)

[32]

28

0.126

2.50

7.22

6.99

10.19

10.05

0.504

0.476

10.52

10.10

0.196

2.00

8.79

8.78

11.51

11.37

0.632

0.602

12.47

11.26

0.256

1.75

9.94

9.79

12.87

12.66

0.727

0.701

13.86

12.47

0.349

1.50

12.11

11.56

15.42

15.23

0.864

0.812

15.93

14.71

0.503

1.25

16.47

16.24

21.47

21.37

1.106

1.081

21.35

20.29

Table 3 Constants in the correlations for calculation of the mean Nusselt numbers, Sherwood numbers and friction factors for the shell side under naturally formed boundary conditions, with staggered arrangement and PT= PL. φ

Nua

PL/do

fa

Sha

C

m

D

n

a

b

0.126

2.50

1.3320

0.4451

2.1168

0.3929

1.6216

-0.25

0.155

2.25

1.3617

0.4562

2.2215

0.3979

1.8357

-0.254

0.196

2.00

1.5146

0.4521

2.3452

0.4038

2.1066

-0.258

0.256

1.75

1.6371

0.4581

2.4831

0.4099

2.5186

-0.266

0.349

1.50

1.9145

0.4569

3.1208

0.3915

3.1984

-0.282

0.503

1.25

2.3844

0.4697

3.7533

0.4162

4.1918

-0.288

29

Table 4 Constants in the correlations for calculation of the mean Nusselt numbers, Sherwood numbers and friction factors for the shell side under naturally formed boundary conditions, with inline arrangement and PT= PL. φ

Nua

PL/do

fa

Sha

C

m

D

n

a

b

0.126

2.50

0.2050

0.6689

0.1373

0.7689

0.5263

-0.192

0.155

2.25

0.2558

0.6502

0.2015

0.7203

0.6352

-0.202

0.196

2.00

0.3449

0.6147

0.3459

0.6473

0.8505

-0.227

0.256

1.75

0.4258

0.5959

0.5887

0.5754

1.1779

-0.257

0.349

1.50

0.5946

0.5597

0.7420

0.5545

2.2086

-0.337

0.503

1.25

1.1762

0.4928

1.5814

0.4732

4.8336

-0.414

30

(a)

H

Air in PL PD

A1

PT

A2

L do y

Solution in

z

W Computational cell

x

(b) Fig.1. Hollow fiber membrane bundle for humidification (MBH): (a) Real photo of the MBH; (b) The schematic.

31

Symmetric boundary

y Air inlet

PL N4 N 1

I

H

A L4

PT/2

E F

Solution

I

Solution

C

B

M4 M1 M2 x

L2 L3

N3

H

G

E

F

Membrane

M3 Membrane

N2

N4

Non-slip boundary D

S

G

Solution

0 L1

N2 N3

Solution

N1

η

A L4 0

B Solution

L1

C

L3

M4

L2

D

Solution

M1

Symmetric boundary

(a)

(b)

Membrane

y*

0.5

0.25

0

0

0.5

1

1.5

2

x* (c) Fig.2. The coordinate system of the periodic computational cell: (a) The physical domain of the computational cell; (b) The rectangular calculating domain of the computational cell; (c) The mesh structure for the physical domain.

32

M3

M2 ε

Start Establish a body-fitted coordinate system; discretize the calculating domain and the governing partial differential equations Solve the continuity, momentum and k-ε Eqs. (1), (18)-(20) and (24)(27); calculate the flow field of the two streams and the k and ε of the air Set the T, ω and X boundary conditions of the two streams Assume the initial temperature (Told*), humidity (ωold*) and mass fraction (Xold*) fields of the two streams Let k=1, z(k)=Δz Solve the energy Eqs. (2) and (21); calculate the temperature fields of the two streams

Update T *= T * old

new

Combined with heat balance Eq. (57) to calculate the membrane temperature, and then calculate a new temperature (Tnew*) fields for both streams No

|Tnew*-Told*|<10-4 Yes Assume an initial moisture emission rate mv through the membrane

Update Told*= Tnew* ωold*=ωnew* Xold*= Xnew*

Solve the mass Eqs. (3) and (22); calculate the concentration fields of the two streams Update mv= mv´

Update k=k+1

Combined with Eq. (58) to calculate a new moisture emission rate mv´ and then calculate new concentration field (ωnew*, Xnew*)for both streams No

|mv´-mv|<10-8 Yes |Tnew*-Told*|<10-4 |ωnew*-ωold*|<10-4 |Xnew*-Xold*|<10-4

No

Yes Solve the Eq. (40); calculate Tb,a*, ωb,a* of the air stream on x-y planes z(k)=kΔz=L

No

Yes Solve the Eq. (41); calculate Tb,a*, ωb,a* of the air stream over the entire fiber length Output velocity, temperature, concentration and turbulent kinetic energy and its dissipation rate fields, and output the f, Nu and Sh End

Fig.3. Flow chart of the solution procedure for the mathematical model. 33

Cold water Dehumidifier

Tso

T

MBH

Screen

Screen

TD,ao ωD,ao

Tao ωao

T H

T H

Screen

Tai ωai

Tsi

PC

T H

Fan

Fresh water tank

DTS

T

Rotameter

T H

Vt

Electric heater

Pump Heat storage tank

Temperature sensor Humidity sensor Hot saline water Cold water Fresh water Data Transmission

Fig.4. Schematic diagram of the membrane-based air humidification-dehumidification desalination system (MHDD).

34

35

Laminar STD k-ε Low Re k-ε Test

30

Nua

25 20 15 10 5

0

100

200

300

400

500

600

Rea Fig.5. Comparison of the mean Nusselt numbers for the air stream by different models and tested values, PL/do= PT/do=1.75, φ=0.256.

35

35 Laminar STD k-ε Low Re k-ε Test

30

Sha

25 20 15 10 0

100

200

300

400

500

600

Rea Fig.6. Comparison of the mean Sherwood numbers for air stream by different models and tested values, PL/do= PT/do=1.75, φ=0.256.

36

1.1 Laminar STD k-ε Low Re k-ε Test

1.0 0.9

fa

0.8 0.7 0.6 0.5 0.4 0

100

200

300

400

500

600

Rea Fig.7. Comparison of the friction factors for the air stream calculated by different models and tested values, PL/do= PT/do=1.75, φ=0.256.

37

y*

0.5

0.25

0

0

0.5

1

1.5

2

x* (a)

0.5 30 40

20 50

40

70

90

0

30 80 70

50 70

0

0.5

20

80

90 80 60 507

60

40

0

70

60

y*

40

80

0.25

8 0 90

20

50 260 50

20

50 30

80

1

1.5

2

x* (b) 7

0.5

9

8

7

7

5 65

10

4

9

7

4 10

6 8

4

y*

10 9

10 11

8

11

87

0.25

9

4

6 9

9

5

0

0

0.5

1

10 9 65 8

8 12

7

10 11

1.5

2

x* (c)

Fig.8. Contours of velocity vectors and dimensionless turbulent kinetic energy (k*) for the air streams: (a) velocity vectors for the air; (b) dimensionless turbulent kinetic energy by STD k-ε model; (c) Dimensionless turbulent kinetic energy by Low Re k-ε model. PT/do= PL/do=1.5, Rea=400.

38

4000 0

0

00

00 1 00

0

0

40

00

3000

0.5

1

0 0 40

00

200

20

5 00

3 00 0

2000 300 0 20

0

20

50

y*

50

0.25

00 0 40 2 0 0

1000

1000

0

0.5

1.5

2

x* (a) 160

0

0

0

12 0

200 160

16

0

16

12

80

0

20

160

22

2000 22

0

80

80

0.25

16 0

22 0

12

y*

22 0 1 20

16

200

80

80 120

0

0.5

0 200

0

0.5

1

1.5

2

x* (b)

Fig.9. Comparison of the contours of dimensionless dissipation rate (ε*) for the air streams: (a) by STD k-ε model; (b) by Low Re k-ε model. PT/do= PL/do=1.5, Rea=400.

39

φ=0.126 φ=0.155 φ=0.196 φ=0.256 φ=0.349 φ=0.503

45 40 35

Nua

30 25 20 15 10 0

100

200

300

400

500

600

Rea Fig.10. Variations of mean Nusselt numbers ( Nua ) for the air stream with Rea at different module packing fractions.

40

φ=0.126 φ=0.155 φ=0.196 φ=0.256 φ=0.349 φ=0.503

45 40 35

Sha

30 25 20 15 10 5

0

100

200

300

400

500

600

Rea Fig.11. Variation of the mean Sherwood numbers ( Sha ) for the air stream with Rea at different module packing fractions.

41

1.6

φ=0.126 φ=0.155 φ=0.196 φ=0.256 φ=0.349 φ=0.503

1.4 1.2

fa

1.0 0.8 0.6 0.4 0.2

0

100

200

300

400

500

600

Rea Fig.12. Variation of friction factors (fa) for the air stream with Rea at different module packing fractions.

42

PT=PL=2.0do PT=2.0do,PL=1.5do PT=1.5do,PL=2.0do PT=PL=1.5do

32 28

Nua

24 20 16 12 8 0

100

200

300

400

500

600

Rea Fig.13. Effects of Rea on the mean Nusselt numbers with different transverse (PT) and longitudinal pitches (PL) of the fibers.

43

32

PT=PL=2.0do PT=2.0do,PL=1.5do PT=1.5do,PL=2.0do PT=PL=1.5do

28

Sha

24 20 16 12 8

0

100

200

300

400

500

600

Rea Fig.14. Effects of Rea on the mean Sherwood numbers with different transverse (PT) and longitudinal pitches (PL) of the fibers.

44

PT=PL=2.0do PT=2.0do,PL=1.5do PT=1.5do,PL=2.0do PT=PL=1.5do

1.2

1.0

fa

0.8

0.6

0.4 0

100

200

300

400

500

600

Rea Fig.15. Effects of Rea on the friction factor with different transverse (PT) and longitudinal pitches (PL) of the fibers.

2

y*

0.25

3 6.5

5

6

6.534 7.5 21 5

4

2 6 7.5 5

7

6.5

1

0.5

4

3

7

6 6.5

5

6

5

6.5

0

6.5

0

0.5

1

1.5

2

x* (a)

6

5

6 41 5 2

0

7

3 42 85

1 7

7

7

6

6

1

8

8

0

6

8

734 25

y*

0.25

7

4

8

7 3 6

5

8

1

6

1 3 42 5

32

0.5

0.5

1

x* 45

1.5

2

(b)

0.5

7. 4 5

0.5

8 7.5

7.5 9 10

1

748.5

5

6.59

6.5

5

1.5

10

0

9

8 190

8

9

10

7.56.5

0

6.5

0.25

7.5

4

y*

8 459

109

58

11 10

6.5

8

7.5 9 10

2

x* (c)

Fig.16. Contours of dimensionless turbulent kinetic energy (k*) of the air streams for different transverse (PT) and longitudinal pitches (PL): (a) PT/do= PL/do=2.0; (b) PT/do=2.0, PL/do=1.5; (c) PT/do=1.5, PL/do=2.0. Rea=400.

Highlights



Conjugate heat and mass transfer in a cross-flow hollow fiber membrane bundle for air humidification (MBH) are investigated considering air side turbulence.



Two turbulence models namely the standard k-ε model (STD k-ε) and the low-Reynolds-number k-ε model (Low Re k-ε) are adopted.



The transverse pitch of the fibers has a more significant impact on the turbulent kinetic energy distribution compared to the longitudinal pitch.



Large values of k and ε are found in the channel center and in the flow separation regions where velocities are higher.

46