nanofiltration plate-and-frame membrane module

Journal of Membrane Science 325 (2008) 580–591

Contents lists available at ScienceDirect

Journal of Membrane Science journal homepage: www.elsevier.com/locate/memsci

Concentration polarization in a reverse osmosis/nanofiltration plate-and-frame membrane module Ana I. Cavaco Morão, Ana M. Brites Alves, Vítor Geraldes ∗ Instituto Superior Técnico, Technical University of Lisbon, Department of Chemical and Biological Engineering, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal

a r t i c l e

i n f o

Article history: Received 8 January 2008 Received in revised form 20 June 2008 Accepted 10 August 2008 Available online 26 August 2008 Keywords: Nanofiltration Reverse osmosis Mass-transfer correlations Computational fluid dynamics Concentration polarization

a b s t r a c t The flow structure and solute concentration distribution in a nanofiltration/reverse osmosis plate-andframe module with radial thin feed channels that have considerable entrance and outlet effects was determined by computational fluid dynamics (CFD). Simulations were performed for binary aqueous solutions, Reynolds (Re) numbers in the range of 64–570 (based on the channel height) and Schmidt (Sc) numbers between 450 and 8900. The CFD simulations showed that both the velocity field and the solute concentration distribution exhibit important 3D effects and, at rather low Reynolds numbers (Re ≥ 118), flow instabilities start to appear in the entrance/outlet regions. However, those instabilities do not affect significantly the average concentration polarization on the membranes surface up to the maximum Reynolds number simulated (Re = 570). The friction factors predicted by CFD were in agreement with the corresponding experimental values for the range of Re numbers investigated. The simulations allowed the determination of a mass-transfer correlation at vanishing mass-transfer rates and a correlation for mass-transfer correction factor. The obtained mass-transfer correlation at vanishing mass-transfer rates was compared with the ones available in the literature as well as with the Sherwood numbers determined by the velocity variation method, using diluted aqueous solutions of glycerol. It was also found that a generalized mass-transfer correction factor correlation for high masstransfer rates, previously developed for membrane modules with 2D configurations, is still valid to predict the average concentration polarization in the module investigated. © 2008 Elsevier B.V. All rights reserved.

1. Introduction The prediction of concentration polarization is required for design and operation of pressure-driven membrane systems such as nanofiltration (NF) and reverse osmosis (RO). As well known, the phenomena of concentration polarization is responsible for the flux decline observed in many applications, due to osmotic pressure raise, increase of resistance to permeation (e.g. gel formation and scaling) and fouling susceptibility. It may also change the selective characteristics of a membrane, for instance due to surface charge variations [1]. Usually, the concentration polarization for binary aqueous solutions is assessed through mass-transfer models which require the knowledge of a mass-transfer coefficient. The mass-transfer coefficient for simple rectangular feed channels that exhibit 2D distribution of velocity and solute concentration, in laminar flow regime, may be determined for instance by the Lêveque corre-

∗ Corresponding author. Tel.: +351 218417511. E-mail address: [email protected] (V. Geraldes). 0376-7388/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.memsci.2008.08.030

lation [2–4], the Grober correlation [3–5] and the Rautenbach correlation [6,7]. These correlations, derived for systems with impermeable walls (i.e. vanishing mass-transfer rates), have been used together with the film theory to predict concentration polarization in plate-and-frame membrane modules, despite the fact that these modules do not have always perfect rectangular feed channels and the flow perturbances in the entrance and outlet regions are nonnegligible. In the absence of adequate masstransfer correlations, indirect experimental methods such as the velocity variation method (VVM) are commonly used to obtain the mass-transfer coefficients, despite their inherent limitations [4,8–10]. Recently, Geraldes and Afonso [2] developed a new approach to predict concentration polarization. With this methodology the conventional mass-transfer coefficient, k␳ – for low mass-transfer rates and impermeable walls – is corrected through a generalized masstransfer correction factor that depends only on the ratio Jv /k␳ . It was found that this generalized correlation for the mass-transfer correction factor is valid for several geometries with 2D velocity fields and concentration distributions both in laminar and turbulent regimes. Nevertheless, the flow in spiral wound and plate-and-frame membrane modules can exhibit a far more complex 3D structure due

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

581

Fig. 1. Geometry of the feed spacer plate (A) and the 3D fluid domain used for the simulations (B and C).

to the insertion of feed spacers and non-negligible entrance and outlet flow effects. For these systems, with 3D velocity fields and concentration distribution, the suitability of the above-mentioned mass-transfer correction factor correlation must be investigated. In the present work, computational fluid dynamics (CFD), which is recognized as a valuable tool for the comprehension and prediction of both hydrodynamics and concentration polarization, e.g. [11,12], was used to predict concentration polarization distribution in a membrane module with a plate-and-frame geometry (LabStak M20 from Alfa Laval; Denmark, Nakskov). This plate-and-frame module was chosen as a case study, because it exhibits 3D velocity fields and significant entrance/outlet effects and, furthermore, is a versatile membrane module widely used in research, e.g. [7,13–19]. Both the Grober [5,13] and the Rautenbach [7] correlations have been applied to predict the mass-transfer coefficient in this module. Nevertheless, the suitability of these equations to predict k␳ for such a complex system with a 3D flow structure was not demonstrated. Thus, for this system, the choice of the mass-transfer correlation remains largely arbitrary and needs to be clarified. In the current work the adequacy of the mass-transfer correlations existing in the literature as well as the suitability of the experimental VVM were evaluated.

2. Case study The cross-flow filtration unit studied is the LabStak M20 plateand-frame module from Alfa Laval (Denmark, Nakskov) suitable for microfiltration, ultrafiltration, NF and RO, with a variable area from 0.036 to 0.72 m2 , composed of a stack of membranes, feed spacer plates (Fig. 1A) and permeate spacer plates which are connected and compressed between two flanges. The two end flanges provide in/outlet connections for cross-flow, i.e. for feed and retentate, and the compression force to put the membrane/plates stack together. Alternating each feed spacer plate, a permeate spacer/carrier plate is mounted, which has flow channels for collecting and conducting permeate to the outlet collector (Fig. 2). The permeate spacer plates are mounted with a membrane sheet on each side, with its active layer outwards, facing the flow channels in the feed spacer plate. The purpose of the feed spacer plates is to create cross-flow feed channel patterns, and to establish a connection in series of the single membrane sheets [20]. Each feed spacer plate has 30 channels with radial flow. In order to reduce the computational requirements, the simulation of the flow and mass-transfer was carried out only in a domain corresponding to 1/30 of the circular plate (shaded region in Fig. 1A).

Fig. 2. Vertical cross-section of the fluid domain evidencing the fluid circulation direction.

582

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

The geometry is depicted in Fig. 1B and C with dimensions matching those of the feed spacer plate. Three distinct regions in the fluid domain were considered, as illustrated in Fig. 1B: inlet/outlet – regions where the fluid is distributed to the several channels of each plate and the retentate is collected; two channels with permeation – regions facing the membrane active layer where permeation takes place; central cylindric region – fluid region that connects the two permeation channels of the fluid domain. 3. Theory and CFD numerical method 3.1. Governing equations and boundary conditions The governing equations are the mass continuity and the momentum equations. For an unsteady flow of an isothermal incompressible Newtonian fluid these equations are given by [21]:

∇ ·U=0





(1)



∂U + ∇ · (U ⊗ U) ∂t

= −∇ P + ∇ · (∇ U + (∇ U)T ) + Sp

(2)

where U is the fluid velocity vector, P is fluid pressure, Sp is the sum of the forces that act on the fluid per unit of volume (excluding the pressure forces) and  and  are the density and viscosity of the fluid, respectively. For an incompressible binary mixture, the solute mass fraction is governed by the solute continuity equation given by ∂ωA + ∇ · (UωA ) = ∇ · (DAB ∇ ωA ) ∂t

membranes the solvent permeates, reducing both the feed flow rate and increasing the average concentration of solute that reaches the following pair of membranes (vide Fig. 2). In order to apply the periodic boundary condition, which significantly reduces the computational demands, a numerical artifact was introduced to keep the average concentration of solutes constant. With this purpose, a very thin layer with 100 ␮m of thickness was added to the outlet surface S2 with a source of solvent balancing the mass of permeated solvent. The introduction of this region ensures the continuity of the mass fluxes between the inlet/outlet of the flow domain simulated. It is worth remarking that the permeate flow rate in each channel was always less than 5% of the feed flow rate. In order to have the fluid circulation along the channel it was necessary to impose a z-momentum source per unit of volume, Spz , in the whole fluid domain that acts as a driving force instead of the pressure difference between the inlet and outlet surfaces S1 and S2. Regarding the boundary condition at the membrane surfaces, it was assumed that there is no coupling between the mass-transfer and the momentum equations. This approach is only completely adequate for membranes with uniform hydraulic permeability and for dilute solutions, with negligible osmotic pressure. Nevertheless, even for the cases where the coupling effect is significant, the average concentration polarization index is only a function of the average permeate flux and does not depend considerably on the local variation of the permeate flux over the membrane as previously demonstrated [2]. For this reason, the mass-transfer correlation developed in this work is expected to be valid also for the cases where the local permeate flux is not uniform.

(3)

where ωA is the solute mass fraction and DAB the mass diffusivity of the solute. The governing equations were solved, with boundary conditions described below, to determine the velocity, pressure and concentration distributions using the CFD code Ansys CFX-10.0 (Canonsburg, USA), which uses the conservative Control Volume Method to discretize the transport equations. The viscosity and the density of the solution as well as the diffusivity of the solutes in water were assumed constant. As explained in Section 2, the transport equations were solved in the domain depicted in Fig. 1B and C, with the following boundary conditions (the names of the different surfaces are displayed in Fig. 1): (1) A translational periodic interface was imposed between the inlet surface S1 and the outlet surface S2. This condition means that the distributions of velocity, concentration and pressure on the surface S1 are equal to the ones of the surface S2. (2) A symmetry boundary condition was applied to the lateral surfaces S3 and S4. This boundary condition imposes that the gradients of all the variables in the direction normal to these surfaces are zero. (3) The normal component of the velocity vector normal to the membrane surfaces S5 and S6 is equal to the permeate flux (Jv ) and the solute flux is these surfaces is zero (the intrinsic rejection coefficient is equal to one). (4) The non-slip condition was imposed on the remaining solid surfaces. The translational periodic interface (boundary condition 1) was implemented because the number of membrane pairs assembled in the plate-and-frame module is variable. If more that one pair of membranes is mounted, in most part of the module the velocity, pressure and concentration distributions exhibit a quasi-periodic behavior. The system does not exhibit a perfect periodicity relatively to the pairs of membranes in series because at each pair of

3.2. Mesh generation The Ansys Design Modeler 10.0 was used to create the geometry and CFX-Mesh 10.0 was used to generate the mesh. The meshing strategy was set to Delaunay surface meshing and advancing front for volume meshing. For the simulation of the flow and mass-transfer, the extremely thin concentration boundary layer requires a large number of control volumes near the membrane surface to solve the concentrations profiles accurately. Thus, a non-uniform mesh with higher density in the thin permeation channels was generated. In order to assess the grid independence of the model output, four distinct meshes were compared with the characteristic depicted in Table 1. More than 70% of the mesh elements were accommodated in the channels contacting the membranes. In the inlet/outlet and central cylindrical regions, the mesh was not so refined, because in these regions solute concentration gradients are low. Therefore, in these regions the mesh was not stretched, the inflated boundaries were not inserted, the maximum body spacing was set to 0.4 mm and the maximum and minimum face spacing edge lengths were 0.2 mm and 0.2 ␮m, respectively. 3.3. Simulation conditions Simulations were performed with CFX-solver 10.0 considering aqueous solutions at 25 ◦ C ( = 8.9 × 10−4 Pa s,  = 997 kg m−3 ) and solute mass fractions ωA0 = 0.001. Six diffusivities were chosen, corresponding to Schmidt numbers in the range 450–8900, covering the characteristic diffusivities of the solutes usually separated by NF and RO. The average cross-flow velocity, v0 , based on the average width of the permeation channel (with an average cross section area of 5.38 mm2 ), was varied by changing the value of Spz between 50,000 and 80,000 kg m−2 s−2 . Simulations were performed for six values of the Reynolds number, based on the height of the permeation channel (h = 0.5 mm) and the average cross-flow velocity, up to 570, considering the flow

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

583

Table 1 Characteristics of the investigated meshes in the membrane channels sub-domaina Mesh 1

Mesh 2

Mesh 3

Mesh 4

Maximum body spacing (mm)

0.4

0.4

0.2

0.2

Face spacing edge lengths Maximum (mm) Minimum (␮m)

0.2 0.2

0.2 0.2

0.1 0. 2

0.1 0.2

Inflated boundary near the membranes surface: Number of layers First layer thickness (␮m) Expansion factor

35 0.2 1.23

100 0.2 1.106

35 0.2 1.23

100 0.2 1.106

Total number elements

1.5 × 106

1.9 × 106

3.3 × 106

5.6 × 106

a

All the meshes were stretched by the factors 0.55, 0.55, 5 in the directions x, y, z, respectively.

rates range of the typical pumps used with this system. The permeation fluxes imposed, for each Re and Sc numbers, were 0.0001, 0.001, 0.002, 0.005, 0.01, 0.02 and 0.05 kg m−2 s−1 , falling within the typical range of fluxes of NF/RO processes. The CFX-solver was first run in steady-state mode with laminar flow. The solver started from the initial conditions of the velocity and pressure set to zero and the solute mass fraction set to 0.001. Using a time step computed by the CFX default time scale algorithm, the solver makes a series of iterations until a given convergence criterion is achieved. The convergence criterion was to set the root mean square of the normalized values of the discretized equations residuals (RMS residuals) equal to 10−4 which is often sufficient for most engineering applications [21]. In addition, the mean value of the concentration polarization was monitored and the simulations were stopped only when steady values were reached. For Re ≥ 120, the solver did not converge in laminar steady-state mode. In order to understand the lack of convergence unsteady simulations were carried out with a time step of 0.25 ms, over a time interval of 0.4 s.

(i.e. without suction), k␳ , through k␳• = k␳ [22]. The mass-transfer correction factor, , depends on the ratio  = Jv /k␳ . Considering the boundary condition of solvent permeation on the membrane, in steady-state and for a intrinsic rejection coefficient equal to one, ¯jAm is given by [2]: ¯jAm = −Jv ¯ Am

(6)

Combining Eq. (4) with the Eq. (6), it follows that Jv ¯ Am − Ao = • ¯ Am k␳

Replacing k␳• in Eq. (7) by the mass-transfer correction factor for high mass-transfer rates, , and inserting the definition of areaaveraged concentration polarization index, ¯ = (¯ Am − A0 )/A0 , the following equation is obtained: ¯ =

1 (/) − 1

using a rearranged version of Eq. (8):

The computed velocity fields and concentrations distributions were subsequently processed using the CFX-Post 10.0 module of Ansys-CFX, in order to determine the average mass-transfer coefficients, the concentration polarization distribution in the membranes and the equivalent pressure drop between the inlet and the outlet of the fluid domain.

k␳ = Jv

¯jAm = −k• (¯ Am − A0 ) ␳

(4)

where ¯jAm is the area-averaged permeation back-diffusion flux of solute A in the direction of the permeation flux vector, which is normal to the membrane surface; A0 is the feed solution concentration and ¯ Am is the area-averaged solute concentration at the membrane surface:



¯ Am =

Am

Am dA Am 

(5)

The black dot “• on the mass-transfer coefficient, k␳• , indicates that the mass-transfer coefficient itself depends on the permeate flux, due to the distortion of the velocity and concentration profiles by the flow through the interface. For permeable walls, with high mass-transfer rates, the value of k␳• is related to the conventional mass-transfer coefficient at vanishing mass-transfer rates

(8)

At vanishing permeation rates, k␳• = k␳ . Thus, under this condition, the mass-transfer coefficient, k␳, can be calculated with ¯ and Jv ,

3.4. Post-processing

3.4.1. Average mass-transfer coefficients The average mass-transfer coefficient, k␳· , is defined by [22]:

(7)

1 ¯



+1

(9)

3.4.2. Calculation of the equivalent pressure drop In the CFD numerical method, resulting from the application of a periodic boundary condition to the inlet and outlet of the fluid domain, the average pressure in both surfaces is the same. As mentioned, the driving force for motion in the numerical method is a momentum source per unit of volume, Spz , in the whole fluid domain, while in the real system, the flow is pressure driven. Nevertheless, at the same flow rate, in both systems the drag force in the z-direction is the same and a balance to the z-component of the forces acting on the fluid leads to (Pin − Pout )eqv Ain = Spz V

(10)

where Pin and Pout are, respectively, the pressure at the inlet and the outlet of the fluid domain, Ain is the area of the inlet surface (73 mm2 ) and V is the volume of the fluid domain (1214 mm3 ). Rearranging Eq. (10), the pressure drop over one pair of membranes, is obtained: (Pin − Pout )eqv . It is assumed that the pressure drop for n pairs of membranes is n × (Pin − Pout )eqv . Regarding the way the whole system is assembled [20], both on the bottom and on the top of the membranes stack, the feed spacer plate contacts directly with the two flanges. Therefore, when the feed solution enters and exits the system it flows through an additional series of channels, shaped like the channels contacting the membranes, i.e. to calculate the pressure drop from the CFD simulations, if n pairs of mem-

584

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

branes are mounted, the pressure drop measured experimentally is equivalent to: (n + 1) × (Pin − Pout )eqv . 3.5. Mass-transfer correlation The values of k␳ obtained for different Re and Sc numbers were used to obtain average mass-transfer correlations at vanishing mass-transfer rates: Sh = a Reb Sc c

(11)

where a, b and c are fitted parameters. Analyzing the physical consistency of Eq. (8), Geraldes and Afonso [2] postulated a new function for (), which is used to obtain a mass-transfer correction factor for high mass-transfer rates:  =  + (1 + c2 c3 )

−c1

(12)

with positive fitted parameters c1 , c2 and c3 . Inserting Eq. (12) into Eq. (8), the following semi-empirical correlation for ¯ is obtained: ¯ =

 −c (1 + c2 c3 ) 1

(13)

The parameters of the previous correlations were obtained through the fitting of Eqs. (11) and (13) to the CFD predictions, using the Nelder-Mead Downhill-Simplex algorithm [23]. 4. Experimental 4.1. Apparatus and materials The NF experiments were performed with NF membranes Desal DK from GE Osmonics (Trevose, PA, USA). Three membrane pairs with an effective membrane area of 0.108 m2 were mounted in the plate-and-frame LabStak M20 unit (Alfa Laval; Denmark, Nakskov). The membranes had an average thickness of 151 ± 6 ␮m and, after the experiments, a thickness of 108 ± 5 ␮m in the regions where they contact the feed spacer plates, meaning that the flat channels were partially hindered. Therefore, the effective channel height was 0.46 mm (assuming that the membrane deformation was irreversible). This value of the channel height was used in all the calculations with experimental data. The thicknesses of the membranes were determined using a digital micrometer with an accuracy of ±1 ␮m. Nanofiltration experiments were performed with an aqueous solution of glycerol with a concentration of 1.25 g/l. This aqueous solution was prepared with deionised water and glycerol, pro analysis grade (supplied by Merck). The feed solution was pumped by a Wanner Engineering (Minneapolis, USA) positive displacement pump (Hydra-Cell Model G/M13) and the temperature was controlled at 15 ± 0.5 ◦ C. The operating transmembrane pressure was measured with an accuracy of ±100 kPa. Concentrate and permeate flow rates were measured by weighing (±0.01 g). The pressure drop was measured with an accuracy of 5 mbar using a multifunction process calibration device from Fluke, model 725. 4.2. Methods 4.2.1. Pressure drop determination The compression of the membranes/plates at 43 MPa ensured a constant channel height for all experiments. The pressure difference between the entrance and the outlet of the plate-and-frame module was measured at nine cross-flow velocities, with three pairs of membranes mounted in the system, using deionised water at 24 ± 1 ◦ C. The pressure sensors were located outside the membrane

stack, in the sockets where the manometers of the system are usually fastened [20]. As the two pressure sensors are not placed in the same horizontal plane, the pressure difference measured was also caused by the gravity force, which acts in the opposite direction of the fluid flow. Thus, the pressure difference due to the gravity force, Pg = gH (H = 6 cm, in the system with three membrane pairs), was subtracted from the values experimentally measured. This pressure drop, caused by the drag force due to the interaction of the fluid with the channels walls, was then used to calculate the friction factor, based on the channel length (L = 51 mm): f =

Pf h v2  L

(14)

0

The friction factor was calculated both using the predicted pressure drop by the CFD simulations and the experimental data, to evaluate the accuracy of the numerical results. 4.2.2. Determination of k␳ by the velocity variation method Preceding first use, the membranes were cleaned in order to remove the preservatives. This cleaning procedure included three steps: (1) rinsing with deionised water; (2) recirculation of deionised water, for 45 min, at 45 ◦ C and pH 10.5, adjusted with NaOH and (3) rinsing with deionised water to neutral pH. After cleaning, the membranes were compacted for 2 h at 4.0 MPa and the pure water flux of the membranes was measured. The hydraulic permeability of the membrane obtained at 15 ◦ C was Lp = 8.5 × 10−12 m s−1 Pa−1 . Total recirculation assays using glycerol solutions were carried out, i.e. the retentate and permeate were both recycled back to the feed tank (8 L) in order to keep the feed concentration constant. Only a minimum volume of permeate (ca. 8 mL), needed for analyses and flux determination, was taken after 10 min of permeation to ensure the system reached a steady state. In order to determine k␳ by the velocity variation method, e.g. [8,9], the observed rejection coefficient of the glycerol, Robs , was measured, at a constant transmembrane pressure, for cross-flow velocities between 0.25 and 1 m s−1 . This procedure was performed for transmembrane pressures of 1.8, 3.0 and 4.0 MPa. The recovery rate, defined by the ratio of the permeate flow rate to the feed flow rate, was always very low and did not exceed 5.8% in any case. The diffusivity of glycerol in water at 15 ◦ C, used to compute the Sh number, is 7.18 × 10−10 m2 s−1 [24]. The feed and permeate concentrations (A0 and Ap ) were determined with a Gilson 133 refractive index detector (Gilson Medical Electronics Inc., Middleton, WI, USA). 5. Results and discussion 5.1. Grid independence study and friction factors The grid independence study was carried out for typical NF conditions (Re = 360 and Jv = 2 × 10−5 m s−1 ) by comparing the areaaveraged concentration polarization index computed with the four meshes having the characteristics depicted in Table 1. The results of the simulations with increasing number of elements (Fig. 3) show clearly that, for all the range of Sc numbers investigated, it is not necessary to refine the mesh above the resolution of mesh 3 in order to predict accurately the average concentration polarization. For this reason, mesh 3, was used in remaining CFD computations. The reliability of computer program to simulate adequately the laminar flow and mass transport, in a rectangular cross section channel with permeable walls, was verified through the comparison of the CFD predictions of concentration polarization (results not shown) with the analytical solution of Sherwood et al. [25].

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

Fig. 3. Concentration polarization index predicted using the numerical method with the four meshes with the characteristics depicted in Table 1.

The accuracy of the CFD numerical method to compute the velocity and pressure fields was evaluated by comparing the friction factors (defined by Eq. (14)) computed by CFD with the experimental values. As shown in Fig. 4, the CFD predictions agree with the experimental results, indicating that the model assumptions, namely the periodic boundary conditions at the inlet/outlet of the feed channels, hold for the range of Reynolds numbers studied. 5.2. Hydrodynamics, solute concentration distribution and concentration polarization At a Re number of 64 the flow is laminar, steady and convergence is achieved in all fluid domain. However, for higher Re numbers the residuals of the momentum and continuity equations increase and convergence (RMS residuals < 10−4 ) is not achieved in steady state modeling due to the transition to unstable laminar flow. The CFD predictions in unsteady-state simulation mode showed that flow whirls arise both in the inlet/outlet regions and in the central cylindrical region. In these regions the velocity and the solute concentrations are fluctuating around a mean value. The time–root mean square (RMS) and time standard deviation (SDV) of velocity, shear stress and concentration were calculated. The ratio between

585

Fig. 4. Variation of the friction factor along the channel with Re.

the time RMS and SDV of these variables was analyzed on the membranes region in order to evaluate the magnitude of the time fluctuations, at Re = 360 and Jv = 2 × 10−5 m s−1 (typical NF operating conditions). These results are depicted in Figs. 5–7 showing that the areas with higher fluctuations of shear stress and concentration polarization over time are coincident as expected. In order to understand the intensity of these instabilities, transient simulations were performed for the set of conditions that exhibit more flow instabilities and for extreme Sc numbers: Re = 570, Jv = 5 × 10−5 m s−1 , Sc = 450 and 8900. The root mean square of the residuals of the momentum and concentration in the membrane channels is lower than 2 × 10−4 (Table 2). In addition, Table 3 shows that the steady-state results are in good agreement with the mean values (RMS) obtained at transient regime despite residuals values. Moreover, for these conditions the ratio SDV/RMS of the shear stress on the membranes is 5.6 and 4.5% for the average polarization index (Sc = 8900). Therefore, although the flow is not steady in the overall fluid domain, the steady-state laminar simulations allow an accurate estimation of concentration polarization on the membranes, despite the fact that full convergence is not achieved in the overall fluid domain. The solutions obtained with the laminar steady-state model were used hereafter. It was always confirmed that a sufficient num-

Fig. 5. Ratio of time standard deviation/time root mean square of velocity (%) in the central plane of the membrane channels and in a longitudinal section of the flow domain at Re = 360, and Jv = 2 × 10−5 m s−1 . (For colour images, the reader is referred to the web version of the article.)

586

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

Fig. 6. Ratio of time standard deviation/time root mean square of shear stress (%) on the membranes surface at Re = 360, and Jv = 2 × 10−5 m s−1 . (For colour images, the reader is referred to the web version of the article.) Table 3 ¯ on the membrane walls and Comparison between the area-average shear stress ( ) area-average concentration polarization (¯ ) obtained for steady state and unsteady−5 1 −1 state simulations (Re = 570, Jv = 5 × 10 m s ) Sc

450 8900

Steady state

Unsteady-state simulationsa

¯ (s−1 )

¯

¯ (s−1 ) RMS

¯ (s−1 ) SDV

¯ RMS

¯ SDV

13,577 13,577

0.621 18.118

13,407 13,407

756 756

0.632 18.347

0.018 0.396

a The root mean square (RMS) and standard deviation (SDV) were computed from the instantaneous area-average values during 0.25 s of simulation time.

Fig. 7. Ratio of time standard deviation and time root mean square of concentration polarization at Re = 360, and Jv = 2 × 10−5 m s−1 and Sc = 890. (For colour images, the reader is referred to the web version of the article.)

ber of iterations were performed to reach the steady-state value of ¯ . The main features of the fluid flow pattern for a Reynolds number of 360 are depicted in Fig. 8. The flow is not simply radial due to the geometry of the feed spacer plate. Mainly due to the inlet and outlet effects there are regions showing pronounced 3D flow characteristics. The flow whirls formed in the inlet region (due its larger cross section) are smoothed out when the fluid enters the

thin channel. The sudden section enlargement in central cylindrical region leads to lower velocities and in the area of the permeation channel behind the central cylindrical region the fluid is stagnated. As fluid goes up, through the central cylindrical region, a complex recirculation pattern is observed but, when it enters the thin channel again, the whirls dragged are dissipated and the flow gets steady, until the section enlargement at the outlet distorts the velocity profile. Clearly the velocity profiles are different in the top and bottom membranes (Fig. 8A and B), not only because of the inlet/outlet and central cylindrical region effects, but also due to the channel cross section expansion and contraction. Fig. 9 shows the complex concentration distribution along the membranes surfaces for typical NF operating conditions. Clearly the change of direction of the fluid coupled with the change of

Table 2 Root mean square of the residuals, in the different regions of the fluid domain, for steady-state laminar simulations, Re = 570, Jv = 5 × 10−5 m s−1 Membrane channels

Central cylindrical region

Inlet

Outlet

Overall fluid domain

1.98 × 10−4 2.15 × 10−4 2.20 × 10−4 2.04 × 10−4 1.36 × 10−4 8.30 × 10−5

2.10 × 10−5 1.70 × 10−5 1.40 × 10−5 8.00 × 10−6 3.00 × 10−6 1.00 × 10−6

1.51 × 10−4 1.20 × 10−4 1.00 × 10−4 5.60 × 10−4 1.80 × 10−4 5.00 × 10−4

2.66 × 10−4 2.13 × 10−4 1.76 × 10−4 9.80 × 10−5 3.10 × 10−5 1.00 × 10−5

1.88 × 10−4 1.92 × 10−4 1.91 × 10−4 1.73 × 10−4 1.21 × 10−4 8.10 × 10−5

7.70 × 10−5 5.20 × 10−5 3.00 × 10−5

7.60 × 10−4 5.77 × 10−4 7.51 × 10−4

1.25 × 10−3 7.77 × 10−4 1.00 × 10−3

1.57 × 10−3 3.57 × 10−4 8.86 × 10−4

5.74 × 10−4 2.57 × 10−4 3.88 × 10−4

Bottom

Top

1.64 × 10−4 1.76 × 10−4 1.78 × 10−4 1.69 × 10−4 1.30 × 10−4 9.40 × 10−5 1.10 × 10−4 7.60 × 10−5 4.10 × 10−5

Sc 450 670 890 1780 4450 8900

vx momentum vy momentum vz momentum

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

587

Fig. 8. Instantaneous velocity contour plot at the center of the membrane channels and in a longitudinal section of the fluid domain. (For colour images, the reader is referred to the web version of the article.)

the cross section area in the inlet/outlet of each feed channel promotes mixing, lowering the concentration polarization (Fig. 9A and B) due to the impingement and suction effects. The complex concentration polarization distribution on the membranes surface results from the complex flow patterns that occur in the feed channel. These complex flow patterns can be observed, for instance, in Fig. 8. It is clear from this figure that the fluid velocity has pronounced variations, not only in the z-direction, but also in the x and y directions. Clearly, the areas with higher time fluctuations (Fig. 7) are also those with lower concentration polarization, showing the effect of the velocity fluctuations on the cleaning of the membrane surface. 5.3. Correlations to predict the average concentration polarization index The mass-transfer coefficients for vanishing mass-transfer rates at the membrane wall, i.e. without distortion of the velocity field and the concentration profiles caused by the permeation flux, were calculated. Thus, in order to determine k␳ , Jv was set to a low value of 10−7 m s−1 ( < 0.02). Thirty-six values of k␳ , based on six Schmidt numbers and six Reynolds numbers, were obtained after post-processing of the CFD results, using Eq. (9). The following mass-transfer correlation was obtained by fitting Eq. (11) to the Sh numbers predicted by CFD: Sh = 0.142Re0.46 Sc 0.37

(15)

This correlation has average and maximum relative errors of 1 and 4%, respectively, in the range of 64 < Re < 570 and 450 < Sc < 8900. The Reynolds number exponent in the mass-transfer correlations for feed channel geometries with 2D distribution of velocity and concentration, at laminar flow regime, is either 0.33 or 0.5 depending whether the velocity profiles are fully developed or developing, e.g. [3,4]. As expected, the Re exponent determined by CFD (Eq. (15)) is within that range, because the velocity profile is not

fully developed in the whole channel (vide Fig. 8). Also, accordingly to the typical mass-transfer correlations at laminar flow regime [2–7], the exponent of the Schmidt number is 0.33. However, the value Sc exponent obtained through the data fitting is somewhat higher, 0.37. The deviation of this exponent from 1/3 must be related to the 3D velocity profiles in the feed channel, caused by the sudden expansion/contraction of cross-section area of the fluid domain at the inlet/outlet of the two feed channels. A similar explanation was already given by Schwinge et al. [26] who point out that in spacer-filled channels the Schmidt number exponent of 1/3 cannot be assumed a priori because the boundary layer development is not similar to an empty channel due to recirculation and reattachment regions. Nevertheless, if the data fitting is carried out setting the Schmidt number exponent to 1/3, the average and maximum relative deviation will increase only slightly to 4 and 13%, respectively. In this case, the following mass-transfer correlation was obtained: Sh = 0.185Re0.46 Sc 0.33

(16)

The experimental method of velocity variation (e.g. [8–10]) was also used to determine an average mass-transfer correlation for the plate-and-frame module in study. This method is based on the description of the concentration polarization by the film theory, which is adequate for moderate permeate fluxes. The experiments with variable cross-flow velocity were carried out with a constant permeate flux and  < 1. Introducing the mass-transfer correlation (Eq. (11)) in a rearranged version of the film theory equation, the following equation is obtained [9,10]: ln

1 − R

obs

Robs



= ln

1 − R R

+

1 Jv h aDAB Sc c Reb

(17)

In this method the exponent b must be imposed and was set to 0.46, as predicted by CFD (Eq. (15)). Fig. 10 displays the plot of ln((1 − Robs )/Robs ) vs. Re−0.46 for the glycerol aqueous solution at different pressures. The experimental results plotted in this way are fairly linear, suggesting the applicability of this method. Nevertheless, as shown in Fig. 10, even at the highest average velocity

588

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

Fig. 9. Instantaneous concentration polarization on the membrane surfaces and concentration distribution in the inlet/outlet regions for typical NF operating conditions. (For colour images, the reader is referred to the web version of the article.)

achieved by the pump (1 m s−1 ), Re−0.46 is quite far from zero, which leads to extrapolation errors in the value of the intercept determined. The intrinsic rejections of the solute was determined by the value of the intercept of these straight lines and used to compute k␳ for each Reynolds number.

Fig. 10. Determination of k␳ using the velocity variation method with a glycerol solution. Extrapolation of ln[(Robs − 1)/Robs ] to Re−0.46 = 0.

In order to obtain a mass-transfer correlation from the experimental results, the value of exponent of the Sc number, c, was set to 0.37, again according to the CFD results. The following correlation was obtained (confidence level 95%, 21 experimental values): Sh = (0.117 ± 0.04)Re0.46 Sc 0.37

(18)

In addition to the extrapolation error, the VVM is based on the assumption that intrinsic rejection coefficient is independent of the feed flow rate. This assumption may not be completely correct, because as the feed flow rate varies at constant transmembrane pressure, the intrinsic rejection coefficient is expected to change slightly due to changes in the concentration at the membrane wall. As this effect is difficult to quantify, it is not possible to fully assess the method error. Besides, as mentioned, to apply this method the Re and Sc exponents must be postulated, which might not always be straightforward as it was shown in this work. The inherent errors associated to this indirect method have already been noticed by other authors, e.g. [9,10]. For these reasons, the comparison of the CFD predictions of the mass-transfer coefficient with the experimental values obtained by the VVM must be done with some caution. The values of Sh obtained by CFD, the correlations obtained by the curve fitting of the CFD results (Eqs. (15) and (16)) and the VVM results (Eq. (18)) as well as the Levêque, Grober and Rautenbach correlations are compared in Fig. 11. The mass-transfer correlations existing in the literature strongly underestimate Sh with average

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

and maximum deviations of 21 and 51% to the Lêvéque correlation; or 35 and 68% to the Rautenbach correlation. The deviations obtained with the Grober correlation are even higher meaning that this correlation is definitely inappropriate for this system. The underestimation of the Sherwood numbers predicted by these correlations, relatively to the CFD predictions, is expected because these correlations do not account for the flow instabilities that occur at the inlet and outlet of the feed channels. In the previous section it was shown that these instabilities contribute to diminish the concentration polarization index and, as a consequence, to increase the average Sherwood number. As already explained, the magnitude of the errors associated to the VVM predictions is difficult to quantify and, for this reason, the agreement between the VVM predictions and the Lêveque and Rautenbach correlations, must be regarded with caution. The mass-transfer correction factor for high mass-transfer rates, , was determined through the fitting of the parameters c1 , c2 and c2 in Eq. (12), minimizing the average relative error between ¯ predicted by CFD and the values calculated by Eq. (13), yielding:  =  + (1 + 0.41.26 )

−1.7

589

Fig. 12. CFD predictions (points) of ¯ vs.  calculated with the mass-transfer correlation (Eq. (15)) together with different equations for .

(19)

This correlation, fitted using the Sh values computed by Eq. (15), is valid for the range of  obtained in the numerical simulations, i.e.  < 8.8. The average and maximum relative errors of ¯ are 3 and 13%, respectively. The suitability of the generalized mass-transfer correction factor correlations, derived for geometries with 2D velocity fields and concentration distributions ([2] and Appendix A), for the system studied in this work was also evaluated. The generalized masstransfer correction factor correlation determined by Geraldes and Afonso [2] is adequate to correct mass-transfer coefficients that are determined for impermeable dissolving walls. For the system investigated in this work, the mass-transfer coefficients were predicted by CFD at vanishing mass fluxes, but with a semi-permeable wall. The computations made in [2] were repeated, as described in Appendix A, for the case of mass-transfer coefficients computed with vanishing mass fluxes and with semi-permeable walls. With this approach the new mass-transfer correction factor correlation, Eq. (A2), was determined. This correlation is valid for 300 < Sc < 3000 and  < 20 and it appears that it can also be generalized for the system under study. Indeed, if Eq. (A2) is used together

with Eq. (15) to predict the average concentration polarization, the relative deviations to the CFD predictions do not increase significantly: average deviations of 5% and a maximum deviation of 19% are achieved for the highest Schmidt number. These results show that the generalized Eq. (A2) can be applied also for a system with 3D velocity profiles and concentration distribution even though the most specific correlation is Eq. (19). These results strengthen the evidence that the correction factors for suction can be taken as independent of the geometry and the flow regime. The average concentration polarization predicted by CFD is compared in Fig. 12 with ¯ calculated using the mass-transfer correlation Eq. (15) together with  obtained by Eq. (19), by the film theory and by Eq. (A2). The film theory predicts ¯ with systematic deviations for  > 2 which increases steeply with . According to these results the film theory may be applied only to  < 2.5, yielding errors below 20%. The observed deviation between the CFD predictions of ¯ and the film theory (Fig. 12) is expected according to previous theories [2,26]. In fact the limitations of the film theory arise only at high values of  due to the suction effects, i.e. the increase of mass-transfer coefficient with permeation flux. These effects are not accurately accounted by the film theory. Aiming for the clarification of the applicability of this extensively used model, Zydney [27] presents a mathematical analysis for the stagnant film assumption, showing that for laminar flow over a flat membrane (assuming constant properties) the film theory is rigorous when  < 2.8. Hereafter, in order to obtain average mass-transfer correlations, the need to perform simulations at several mass-transfer rates, i.e. permeation velocities, is eliminated reducing significantly the computing time required. In the present work, it would represent a reduction from 252 to 36 simulations, given that only the results for low permeation velocities at the membrane wall are needed. 6. Conclusions

Fig. 11. Comparison of the Sherwood numbers predicted by CFD (Jv = 10−7 m s−1 ) with Eqs. (15) and (16), and the correlations of Grober 0.5 Sc 0.33 (D /L)0.5 ), 0.33 Sc 0.33 ) Lévêque (ShDh = 1.47ReD and (ShDh = 0.664ReD h h

h

1/3

Rautenbach (ShDh = [3.663 + 1.613 ReDh Sc(Dh /L)] ). The dimensionless numbers of these correlations are defined in the list of symbols.

CFD proved to be a reliable tool to derive mass-transfer correlations for a plate-and-frame module with a 3D flow structure and highly non-uniform concentration polarization. It was observed that at rather low Reynolds numbers (Re ≥ 118) flow instabilities start to appear in some regions of the system. However, it was shown that those time fluctuations do not affect significantly the average concentration polarization on the membranes surface up to the maximum Reynolds number simulated (Re = 570).

590

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

Through the CFD predictions of concentration polarization for laminar flow a mass-transfer correlation for vanishing masstransfer rates was obtained for the plate-and-frame module LabStak M20: Sh = 0.142Re0.46 Sc0.37 . The mass-transfer correlations available in the literature for rectangular thin channels overestimate the average concentration polarization index in this system because the instabilities of the flow in the entrance regions contribute to clean the membranes. Nevertheless, the Lêvéque correlation may be used if an average overestimation of 20% (maximum 51%) in the average concentration polarization index is considered adequate for the range of Re numbers investigated. On the other hand, the Grober correlation cannot be applied to this system because it leads to quite high deviations. The generalized correction factor correlation for high masstransfer rates derived for membrane modules with 2D geometries by Geraldes and Afonso [2] and modified in this work for masstransfer coefficients obtained for semi-permeable walls was found to be appropriate also for module investigated, which has a more complex 3D geometry. Therefore, based on this fact, the number of CFD simulations required to develop new correlations to predict the average concentration polarization in other NF/RO modules may be significantly reduced (from 252 to 36, in the present work). Even though, for aqueous solutions with low Sc numbers and permeation fluxes typical of NF processes, the film theory can be applied to predict the concentration polarization. For processes dealing with higher Sc numbers and higher fluxes a generalized mass-transfer correction factor correlation [2] should be used instead. Acknowledgement The PhD grant SFRH/BD/18496/2004 from the Portuguese Foundation for Science and Technology is gratefully acknowledged by Cavaco Morão, A.I. Appendix A. Generalized correction factor correlation In a previous paper [2], Geraldes and Afonso developed a generalized mass-transfer correction factor correlation for NF and RO, which enables the expedite computation of concentration polarization for binary solutions, based on the Sherwood number, Sh, for vanishing mass-transfer rates in a wide range of operation conditions. In that work, the Levêque correlation was used for obtaining Sh in a rectangular slit with developed laminar flow. As well know, the Levêque mass-transfer correlation was developed – by analogy with heat-transfer – for the case of a boundary condition of constant solute concentration at the channel walls. If instead of that boundary condition, another one is used which assumes unhindered transport of solute and solvent across the membrane with an intrinsic rejection coefficient R, a different mass-transfer correlation for vanishing mass-transfer rates should be obtained. Using the same Fortran code used in Ref. [2], CFD computa −5 tions  were performed for vanishing permeate fluxes, i.e. j = 10 , and for the remaining dimensionless numbers of Table 1 j,max displayed in Ref. [2]. The Sh correlation obtained with these new numerical simulations is Sh = 1.56(Re Sc h/l)

1/3

(A1)

This correlation is very similar to the Levêque correlation, with only a small difference in the first constant factor. Using this Sh correlation, a new mass-transfer correction coefficient was obtained, based on the same approach described in Ref. [2]:  =  + (1 + 0.31.4 )

−1.7

(A2)

These two correlations should be used whenever the computations of Sh are based on the boundary condition of unhindered transport of solute and solvent across the membrane with an intrinsic rejection R.

Nomenclature area of the inlet surface (m2 ) membrane area (m2 ) diffusivity of the solute A (m2 s−1 ) hydraulic diameter; for a narrow slit geometry, Dh ∼ 2 h (m) f friction factor for the fluid flow in the module channel (Eq. (14)) g acceleration of gravity (9.8 m s−2 ) h channel height (m) H height between the two pressure sensors (m) j diffusive flux (kg m−2 s−1 ) permeation flux (m s−1 ) Jv k␳ mass-transfer coefficient (m s−1 ) L channel length (m) hydraulic permeability (m s−1 Pa−1 ) Lp P transmembrane pressure (Pa) P transmembrane pressure (Pa) Pl pressure loss along the channel (Pa) Pg pressure loss due to gravity (Pa) Pin , Pout , average pressure at the inlet and the outlet of the fluid domain, respectively (Pa) R intrinsic rejection coefficient, R = (¯ Am − Ap )/¯ Am Re Reynolds number, Re = v0 h/ Reynolds number based on the hydraulic diameter, ReDh ReDh = v0 Dh / observed rejection, Robs = (A0 − AP )/A0 Robs Sc Schmidt number, Sc = /DAB Sh Sherwood number, Sh = k␳ h/DAB ShDh Sherwood number based on the hydraulic diameter, ShDh = k␳ Dh /DAB Sp sum of the forces that act on the fluid per unit of volume (kg m−1 s−2 ) Spz momentum source per unit of volume in the z direction (kg m−1 s−2 ) t time (s) U fluid velocity vector (m s−1 ) V volume of the fluid domain (m3 ) v0 average cross-flow velocity at the membrane channel (m s−1 ) x x-coordinate (m) y y-coordinate (m) z z-coordinate (m) Ain Am DAB Dh

Greek letters ¯ area-averaged concentration polarization index, ¯ = (¯ AM − A0 )/A0  

 j

 A

¯  ωA

viscosity (kg m−1 s−1 ) mass-transfer coefficient correction factor,  = k␳• /k␳ dimensionless pure water flux (˘j = Lp P/v0 ) density (kg m−3 ) mass concentration of solute A (kg m−3 ) area-averaged shear stress (kg m−1 s−2 ) ratio between Jv and k␳ ,  = Jv /k␳ solute mass fraction

A.I. Cavaco Morão et al. / Journal of Membrane Science 325 (2008) 580–591

Subscripts A solute FT film theory m at the membrane surface P permeate 0 bulk flow in the feed side Superscripts • indicates that the mass-transfer coefficient depends on the mass-transfer rate a, b, c positive fitting parameters in Eq. (11) c1 , c2 , c3 positive fitting parameters in Eq. (12)

References [1] M.D. Afonso, G. Hagmeyer, R. Gimbel, Streaming potential measurements to assess the variation of nanofiltration membranes surface charge with the concentration of salt solutions, Sep. Purif. Technol. 22/23 (2001) 529. [2] V. Geraldes, M.D. Afonso, Generalized mass-transfer correction factor for nanofiltration and reverse osmosis, AIChE J. 52 (10) (2006) 3353. [3] V. Gekas, B. Hallström, Mass transfer in the membrane concentration polarization layer under turbulent cross flow. I. Critical literature review and adaptation of existing correlations to membrane operations, J. Membr. Sci. 30 (1987) 153. [4] A.I. Shäfer, A.G. Fane, T.D. Waite, Nanofiltration—Principles and Applications, Elsevier, Oxford, 2005. [5] R.F. Madsen, Hyperfiltration and Ultrafiltration in Plate and Frame Systems, Elsevier, New York, 1989. [6] R. Rautenbach, R. Albrecht, Membrane Processes, Wiley, NY, 1989. [7] H.C. van der Horst, J.M.K. Timmer, T. Robbertsen, J. Leenders, Use of nanofiltration for concentration and demineralization in the dairy industry: model for mass transport, J. Membr. Sci. 104 (1995) 205. [8] A. Bouchoux, H. Roux-de Balmann, F. Lutin, Nanofiltration of glucose and sodium lactate solutions—variations of retention between single and mixed solute solutions, J. Membr. Sci. 258 (2005) 123. [9] G.B. van der Berg, I.G. Rácz, C.A. Smolders, Mass transfer coefficients in cross flow ultrafiltration, J. Membr. Sci. 47 (1989) 25.

591

[10] M.S.H. Bader, J.N. Veenstra, Analysis of concentration polarization phenomenon in ultrafiltration under turbulent flow conditions, J. Membr. Sci. 114 (1996) 139. [11] R. Ghidossi, D. Veyret, P. Moulin, Computational fluid dynamics applied to membranes: state of the art and opportunities, Chem. Eng. Process. 45 (2006) 437. [12] D.E. Wiley, D.F. Fletcher, Techniques for computational fluid dynamic modelling of flow in membrane channels, J. Membr. Sci. 211 (2003) 127. [13] A.I.C. Morão, A.A.M. Brites, M.C. Costa, J.P. Cardoso, Nanofiltration of a clarified fermentation broth, Chem. Eng. Sci. 61 (2006) 2418. [14] M. Mänttäri, A. Pihlajamäki, M. Nyström, Comparison of the nanofiltration and tight ultrafiltration membranes in the filtration of paper mill process water, Desalination 149 (2002) 131. [15] F. Lipnizki, J. Boelsmand, R.F. Madsen, Concepts of industrial-scale diafiltration systems, Desalination 144 (2002) 179. [16] G. Capar, U. Yetis, L. Yilmaz, Reclamation of printing effluents of a carpet manufacturing industry by membrane processes, J. Membr. Sci. 277 (2006) 120. [17] M. Afonso, Surface charge on loose nanofiltration membranes, Desalination 191 (2006) 262. [18] M.R. Teixeira, M.J. Rosa, The impact of the water background inorganic matrix on the natural organic matter removal by nanofiltration, J. Membr. Sci. 279 (2006) 513. [19] M. Mänttari, K. Viitikko, M. Nyström, Nanofiltration of biologically treated effluents from the pulp and paper industry, J. Membr. Sci. 272 (2006) 152. [20] Operating manual DSS Module Type LabStak® M20-0.72-PSO, DSS file 400016, Denmark, Nakskov, 1999. [21] ANSYS CFX-Solver help navigator, Release 10.0: Modelling, Canonsburg, USA, document built on July 2005. [22] B.R. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, John Wiley & Sons, Madison, WI, 1960. [23] L.I.R. Garcia, Nelder Mead Optimization Algorithm in: Optimization Tool and Equations Solver for EXCEL, version 2.0.1, Foxes Team Eds, 2006. Web site: http://digilander.libero.it/foxes/SoftwareDownload.htm. [24] P.W. Atkins, Physical Chemistry, 6th edn., University Press, Oxford, 1994. [25] T.K. Sherwood, P.L.T. Brian, R.E. Fisher, L. Dresner, Salt concentration at phase boundaries in desalination by reverse osmosis, Ind. Eng. Chem. Fundam. 4 (2) (1965) 113. [26] J. Schwinge, D.E. Wiley, D.F. Fletcher, Simulation of the flow around spacer filaments between channel walls. 2. Mass-transfer enhancement, Ind. Eng. Chem. Res. 41 (2002) 4879. [27] A.L. Zydney, Stagnant film model for concentration polarization in membrane systems, J. Membr. Sci. 130 (1997) 275.