- Email: [email protected]
Simulation of heavy metal extraction in membrane contactors using computational fluid dynamics
Desalination 281 (2011) 422–428
Contents lists available at SciVerse ScienceDirect
Desalination journal homepage: www.elsevier.com/locate/desal
Simulation of heavy metal extraction in membrane contactors using computational fluid dynamics Azam Marjani ⁎, Saeed Shirazian Islamic Azad University, Farahan Branch, Department of Chemical Engineering, Farahan, Iran
a r t i c l e
i n f o
Article history: Received 19 March 2011 Received in revised form 12 August 2011 Accepted 13 August 2011 Available online 9 September 2011 Keywords: Heavy metals Extraction Numerical simulation Computational fluid dynamics
a b s t r a c t A new approach based on momentum and mass transfer for modeling of wastewater treatment was developed in this study. The model is based on solving governing equations for solutes using finite element method (FEM). Wastewater treatment medium considered in this study is a microporous membrane contactor made of polypropylene which is a hydrophobic polymer. Modeling was carried out for separation of water contaminants. With solving conservation equations, velocity and concentration distributions were obtained for the solute in the membrane contactor. Modeling findings showed high removal efficiency for copper extraction. The effect of process parameters on the extraction of copper was investigated. Furthermore, simulation results confirmed that membrane contactors are efficient in the removal of water contaminants. © 2011 Elsevier B.V. All rights reserved.
1. Introduction Nowadays, development of clean technologies for removal of heavy metals from water and wastewaters has been of great importance. From environmental perspectives, it is of vital importance to remove heavy metals from wastewaters to prevent environmental problems. Recently, a non-dispersive extraction process utilizing a porous membrane as contactor has attracted many attentions. In this type of process, porous membrane is used as a contact medium between aqueous and organic phases. Membrane contactor is constituted of two parts, i.e. lumen and shell side. Both aqueous and organic solutions flow continuously, the first one in the lumen of fibers and the second in the shell side. Both phases make interface through pores of the fiber wall [1]. Membrane contactors as extraction devices have become a subject of great interest. In these contactors, the membrane acts as a physical barrier between the two phases, i.e. feed and solvent. Membrane contactors provide high surface/volume ratio. In addition, in membrane contactors, the interfacial mass transfer area is relatively high, whereas the hydrophilic or hydrophobic nature of the membrane determines the position of the interface between the feed and the solvent. A major part of the interest towards membrane contactors is due to their capability in setting a dispersion free contact. In addition, the velocities of both phases can be chosen independently, while neither flooding nor unloading problems arise [2,3].
⁎ Corresponding author. E-mail address: [email protected] (A. Marjani). 0011-9164/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.desal.2011.08.032
Extraction of heavy metals using organophosphorus extractants is an important process in hydrometallurgical industries [4]. Conventional solvent extraction processes are carried out using packed towers, mixer-settlers, etc. which try to maximize the contact area of two phases for separation operation. These contactors have some disadvantages such as formation of stable emulsions, thereby inhibiting the phase separation and product recovery; avoid using liquids having similar densities and foaming. Additional limitations present in packed towers include loading requirements and flooding restrictions [5]. Membrane contactors can overcome the drawbacks of conventional solvent extraction processes. Therefore, extraction of heavy metals can be accomplished using membrane contactors. Some researchers have examined the applicability of these contactors for Meal extraction. Juang et al. [5,6] investigated the extraction of copper (II) using di(2-ethylhexyl)phosphoric acid (D2EHPA) in a flat sheet and hollow-fiber membrane contactor. Their results proved the applicability of membrane contactors for removal of copper (II). They also developed a mass transfer model for prediction of membrane extractor performance. Their model was based on resistancein-series model. Trtic-Petrovic et al. [7] used membrane-based solvent extraction for the separation of thallium (III) from chloride-containing acidic solutions with butyl-acetate. They developed a mass transfer model and showed that the existence of concentration boundary layer (CBL) in the aqueous phase inside the hollow fibers is the main limiting factor in TI(III) transport through the interfacial area. Recently Guo and Ho [8] developed a mass transfer model for simulation of Cu 2+ extraction in a membrane contactor. They solved continuity equation for transport of Cu 2+ across the membrane. They
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
used analytical method to solve the governing equation. Parabolic velocity profile for velocity in the contactor was assumed and axial diffusion was neglected which is not a good approximation especially for low flow rates. There are two approaches for modeling of metals extraction in membrane contactors. The first approach is based on resistances-inseries model. This model considers three mass transfer resistances in series for transport of metal from feed to solvent phase. The mass transfer coefficients for all phases are obtained experimentally using correlations which are not very accurate. The second approach is based on solving the conservation equations including mass, momentum and energy for metal in the three sections of membrane contactor, i.e. lumen, membrane and shell side. For this approach numerical method is applied for solving the equations to predict the performance of membrane extractor. There is a definite need for a comprehensive model that can provide a general simulation of the extraction of heavy metals in membrane contactor [9–11]. The main purpose of the current study is to develop and solve a mathematical model for prediction of metal extraction in membrane contactors. The governing equations are then solved numerically using CFD techniques. The velocity field is determined using Navier– Stokes equations. The aim of this simulation is to predict the performance of membrane contactor for extraction of copper.
2. Theory 2.1. Model equations Fig. 1 shows a tubular membrane contactor for extraction of copper. The feed phase containing aqueous solution of Cu 2+ flows inside the inner tube (lumen) while the extractant solution flows countercurrently inside the shell side. Non-wetted mode is assumed for numerical simulation. It is assumed that the extractant-containing phase penetrates the membrane pores and wets them due to
423
hydrophobicity of the membrane. An interface is established at the pores' mouth adjacent to the lumen side. At the interface, chemical reaction occurs between metal ion (Cu 2+) and organometallic complex is formed and diffuses through the membrane pores and is transferred to the shell side. Therefore, the driving force in this process is concentration difference of solute. The mass transfer model is built considering the following assumptions: (1) (2) (3) (4) (5)
Steady state and isothermal conditions Fully developed parabolic velocity profile inside fiber Laminar flow for two phases in the contactor The membrane is completely extraction-phase filled The chemical reaction is very fast and achieves equilibrium.
The main equation that describes the transfer of solute (Cu 2+) from aqueous phase to organic phase is continuity equation. The differential form of continuity equation for solute may be written as [12]: ∂Ci þ ∇⋅ð−Di ∇Ci þ Ci V Þ ¼ Ri ∂t
ð1Þ
where Ci denotes the concentration of solute (mol/m 3), Di denotes its diffusion coefficient (m 2/s), V the velocity vector (m/s) and Ri denotes the reaction term (mol/m 3.s). The velocity vector can be expressed analytically or obtained by coupling a momentum balance to the equation system. This equation is also called convection and diffusion equation. The diffusion coefficient for the dissolved species accounts exclusively for the interaction between the solute and the solvent. In the Maxwell–Stefan convection–diffusion application mode, the interaction between the different dissolved species is also taken into account. Eq. (1) is the main equation of mass transfer. This equation should be solved numerically to obtain the concentration distribution of solute in the membrane contactor. To solve the continuity equation Eq. (1), velocity distribution is needed. Velocity distribution is obtained by solving the momentum equation, i.e. Navier–Stokes equations. Therefore, the momentum and continuity equation should be coupled and solved simultaneously to obtain the concentration distribution of solute. The Navier–Stokes equations describe flow in viscous fluids through momentum balances for each of the components [12]: ∂V T −∇·η ∇V þ ð∇V Þ þ ρðV:∇ÞV þ ∇p ¼ F ∂t ∇:V ¼ 0
ρ
ð2Þ
where η denotes the dynamic viscosity (kg/m.s), V the velocity vector (m/s), ρ the density of the fluid (kg/m 3), p the pressure (Pa) and F is a body force term (N). The continuity equation Eq. (1) is developed for all subdomains of the membrane contactor, i.e. lumen, membrane and shell sides. It may be written as: " # ∂2 Ci 1 ∂Ci ∂2 Ci ∂C Di þ þ 2 ¼ Vz i−shell : r ∂r ∂z ∂r 2 ∂z
ð3Þ
The velocity distribution in the lumen side is assumed to follow Newtonian laminar flow [12]: 2 r VZ−tube ¼ 2u 1− r1
Fig. 1. Tubular membrane contactor for extraction of Cu2+.
ð4Þ
where u is average velocity in the lumen side and r1 is inner radius of lumen (see Fig. 2).
424
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
2.2. Numerical solution of model equations
Fig. 2. Convergence behavior for numerical simulation.
The velocity in the shell side is characterized using Navier–Stokes equations. Eq. (3) should be scaled because of large difference between membrane radius and length. The equation is scaled in the z-direction. A scaling factor is defined as: δ ¼ z=S
ð5Þ
where S is scale factor. Applying the scale factor results in: " # ∂2 Ci 1 ∂Ci ∂2 Ci ∂C Di þ þ ¼ Vz i−shell : r ∂r S2 ∂δ2 S∂δ ∂r 2
ð6Þ
Arranging Eq. (6) results in:
SDi
∂2 Ci 1 ∂Ci Di ∂2 Ci ∂C þ SDi ¼ Vz i−shell : þ 2 r ∂r S ∂δ2 ∂δ ∂r
ð7Þ
Therefore, anisotropic diffusivity is needed for numerical simulation. The diffusivity in r-direction is SDi while the diffusivity in z-direction is Di S. Boundary conditions for numerical simulation are listed in Table 1. Where m is partition coefficient of solute between aqueous phase and organic phase. C1, C2, and C3 are concentration of metal in the lumen, membrane, and shell side respectively. The boundary condition for the convective flux assumes that the mass passing through this boundary is convection-dominated. In other words, it assumes that the mass flux due to diffusion across this boundary is zero. Diffusion mechanism (V = 0) is assumed for mass transfer of solute within the membrane. The chemical reaction for all sections is assumed to be zero (Ri = 0).
Table 1 Boundary conditions of mass transfer equations. Position
Lumen
Membrane
Shell
z = 0 (Inlet) z = L (Outlet) r = r1 r = r2 r = r3
Ci = C0 Convective flux C1 = C2/m – –
Insulated Insulated C2 = C1*m C2 = C3 –
Convective flux Ci = 0 – C3 = C2 ∂Ci ¼ 0 (symmetry) ∂r
The model equations with the appropriate boundary conditions were solved numerically using COMSOL software. This software employs finite element method (FEM) for numerical solutions of model equations. The finite element analysis is combined with adaptive meshing and error control using numerical solver of UMFPACK [9]. This solver is well suited for solving stiff and non-stiff non-linear boundary value problems [9]. An IBM-PC-Pentium4 (CPU speed is 2800 MHz) was used to solve the set of equations. The computational time for solving the set of equations was about 2 min. Table 2 shows parameters used in the numerical simulation. The maximum number of iterations was set at 25 and the relative tolerance was 1E−6. After about 9 iterations the system reaches convergence. Fig. 2 shows the convergence behavior for numerical simulation. Fig. 2 implies that the system is not stiff and provides confidence in the solution process and ensures consistency. Fig. 3 shows a segment of the mesh used to determine the solute transport behavior in the membrane contactor (MC). It should be pointed out that the COMSOL mesh generator creates triangular meshes that are isotropic in size. A large number of elements are then created with scaling. A scaling factor of 120 has been employed in z-direction due to large difference between r and z. COMSOL automatically scales back the geometry after meshing. This generates an anisotropic mesh around 4397 elements. Adaptive mesh refinement in COMSOL, which generates the best and minimal meshes, was used to mesh the contactor geometry. For numerical simulation, it is important to select a mesh that minimizes the error in the calculations. It is not easy to specify mesh sizes manually to minimize the error. The algorithm must resolve the solution in great detail only on small portions of the domain. Adaptive mesh generation identifies the regions where high resolution is needed and produces an appropriate mesh. The adaptive solver works on linear or nonlinear stationary PDE problems (including parameter-dependent problems) or eigenvalue PDE problems using the coefficient or general solution forms. The adaptive mesh refinement performs the following iterative algorithm for numerical simulation of model's equations: 1. Solve the problem on the existing mesh using the stationary solver. 2. Evaluate the residual of the PDE on all mesh elements. 3. Estimate the error in the solution on all mesh elements. The computed error estimate is really just an error indicator because the estimate involves an unknown constant. 4. Terminate execution if it has made the requested number of refinements or if it has exceeded the maximum number of elements. 5. Refine a subset of the elements based on the sizes of the local error indicators. 6. Repeat from Step 1 [13].
Table 2 Parameters used in the simulations [5,6]. Extractor length (cm) Fiber o.d. (mm) Fiber i.d. (mm) Module i.d. (mm) Scale factor Kex Porosity (ε) Tortuosity (τ) DCu2+-water (m2/s) Dcomplex-organic (m2/s) Dcomplex-membrane (m2/s) Extractant concentration (mol/m3) Inlet concentration of solute (mol/m3)
30 20 22 30 120 4.4E−4 0.67 1.5 3.25E−8 2E−10 Dcomplex-organic (ε/τ) 30 1
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
425
(extractant) flows from the other side (z = L) where the concentration of solute is assumed to be zero. As the feed flows through the lumen side, solute (Cu 2+) is transferred towards the membrane due to the concentration difference. At the membrane-aqueous phase interface chemical reaction occurs and metallic complex is formed. The formed complex diffuses through the membrane pores which is filled with the organic phase and reaches at the shell side. Finally, the complex is swept by the moving extractant and leaves the extractor. Fig. 4 confirms that the concentration decreases along the extractor because of the extraction of solute. Fig. 4 also indicates that increasing feed flow rate decreases the solute removal in the contactor. Increasing feed flow rate decreases residence time of feed phase in the contactor which in turn results in lower mass transfer of solute [14,15]. The extraction efficiency of solute which is defined the ratio of the solute transfer from the feed phase to extraction phase to total solute in the initial feed phase is calculated using numerical simulation. The extraction efficiency for solute is calculated 85.8% for extraction of Cu 2+ using membrane extractor. 3.2. Contours of concentration in the lumen side
Fig. 3. Magnified segment of the mesh used in the numerical simulation (There are 4397 elements in total for the whole MC domain; z-direction scale factor= 120; the three domains from left to right are lumen side, membrane and shell side, respectively).
To investigate the concentration of solute in the lumen side of the membrane contactor, contours of solute concentration is calculated. Fig. 5 shows the contours of solute concentration in the lumen side where the feed phase flows with laminar flow. As it can be seen from Fig. 5, concentration gradient is significant at the regions near
3. Results and discussion 3.1. Concentration profile of solute along the membrane contactor Dimensionless concentration profile (C/C0) of solute is shown in Fig. 4. The concentration profile is obtained in the lumen side of the membrane contactor where the aqueous feed flows. It should be pointed out that the flow pattern in the contactor (extractor) is parallel and counter current. The feed phase including aqueous solution of Cu 2+ flows from one side of the contactor (z = 0) where the concentration of solute is the highest (C0), whereas the organic phase
Fig. 4. Axial concentration distribution of solute along the extractor. Organic flow rate = 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature = 298 K, Pressure = 1 atm.
Fig. 5. Contours of concentration in the lumen side. Feed flow rate = 1E−8 m3/s, organic flow rate = 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature= 298 K, Pressure= 1 atm.
426
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
the inlet of lumen side. High mass transfer flux occurs in these regions of lumen because of higher driving force. As the feed flows along the lumen side, it is transferred to the shell side and causes lower concentration gradient at the regions near the lumen outlet. 3.3. Profile of concentration in the shell side Fig. 6 indicates concentration profile of solute in the shell side of contactor. Obviously, concentration difference at a narrow region near the fiber wall is sharp and significant. This is an evidence for the formation of concentration boundary layer in the shell side. At the regions far from the fiber wall, solute concentration is zero. It reveals that the solute cannot be transferred into the distances far from fiber wall and is swept by the solvent flow. Another reason for such behavior could be attributed to this fact that the solubilizing capacity of the extractant is high and partition coefficient of solute between feed and organic phase is relatively high. 3.4. Streamline of solute concentration in the lumen side Streamline of solute concentration in the lumen side of membrane contactor is illustrated in Fig. 7. It is seen that the solute is transferred from the bulk of feed phase towards the membrane in the r-direction because of concentration difference. The mass transfer mechanism in the r-direction is governed by molecular diffusion. It means that the diffusion is the predominant in the radial direction while in axial direction the convection is predominant. It is clearly seen that at the regions near the lumen center, convection in the z-direction is significant and causes mass transfer of solute towards the lumen outlet. However, to increase the removal rate of solute, concentration gradient in r-direction should be enhanced. This causes the mass transfer of solute to membrane and shell side. This can be accomplished by choosing strong and effective extractant. 3.5. Velocity field in the extractor The contours of velocity in the shell side of membrane contactor are illustrated in Fig. 8. Fig. 9 also indicates the profile of organic phase in the shell side at different axial positions. The velocity field in the shell side of membrane contactor was simulated by solving the Navier–Stokes equations. As it can be seen from Figs. 8 and 9, the velocity profile is almost parabolic with a
Fig. 7. Streamline of solute concentration in the lumen side. Feed flow rate = 1E −8 m3/s, organic flow rate = 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature = 298 K, Pressure = 1 atm.
mean velocity which increases with the membrane length because of continuous fluid permeation. Fig. 9 also reveals that at the inlet regions in the liquid side, the velocity is not developed. After a short distance from the inlet, the velocity profile is fully developed (see Fig. 9). As observed, the model considers the entry effects on the hydrodynamics of fluid flow in the shell side. Furthermore, Fig. 8 reveals the effect of viscous forces at the region adjacent the shell wall. The effect of viscous force in these regions is significant and causes zero velocity (no slip) at the shell wall. 3.6. Effect of porosity-to-tortuosity ratio In membrane contactors, porosity and tortuosity of membrane is of great importance. These parameters have effect on the extraction efficiency of solute [10]. To investigate the effect of porosity and tortuosity on the performance of process, porosity-to-tortuosity ratio varied between 0.05 and 0.80 has been considered in the simulations. Fig. 10 shows the effect of the porosity-to-tortuosity ratio on the extraction efficiency of copper. The removal of copper has grown about 7 times when the latter ratio increases from 0.05 to 0.80. The effective diffusion coefficient (Dcomplex-membrane) of the membrane is calculated through the porosity and tortuosity of the membrane [10]: Dcomplex−membrane ¼ Dcomplex−organic ðε=τÞ:
Fig. 6. Concentration profile of solute in the shell side. Feed flow rate= 1E−8 m3/s, organic flow rate= 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature= 298 K, Pressure = 1 atm.
ð8Þ
As it can be seen from Eq. (8), the effective diffusion coefficient is a function of porosity and tortuosity of membrane. As the porosity/tortuosity ratio increases, the membrane diffusivity and thus the mass transfer of complex through the membrane increases; i.e. the extraction efficiency of solute increases. Indeed, when the porosity/tortuosity ratio increases, the membrane mass transfer resistance decreases. Therefore, the total resistance to the mass transfer of solute becomes lower.
Extraction efficiency of solute (%)
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
427
90 80 70 60 50 40 30 20 10 0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Porosity / tortuosity ratio Fig. 10. Effect of porosity-to-tortuosity ratio on extraction efficiency of solute. Feed flow rate = 1E−8 m3/s, organic flow rate = 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, extractant concentration = 30 mol/m3, Temperature = 298 K, Pressure = 1 atm.
4. Model validation To validate the modeling findings obtained from mass transfer model developed here, the simulation results were compared with the experimental data obtained from literature for extraction of copper using a hollow-fiber membrane contactor. The experimental data are for recycling mode of operation. An additional equation for feed tank is required for simulation of recycling mode. The mass balance equation for feed tank containing aqueous solution of copper (solute) can be written as follows: υ
Fig. 8. Contours of velocity in the shell side. Feed flow rate = 1E−8 m3/s, organic flow rate = 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature = 298 K, Pressure = 1 atm.
Fig. 9. Velocity profile in the shell. Feed flow rate = 1E−8 m3/s, organic flow rate = 1E −7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature = 298 K, Pressure = 1 atm.
dC tank ¼ QC jZ¼L −QC tank dt
@t ¼ 0;
Ctank ¼ Cint
ð9Þ ð10Þ
where Qis volumetric flow rate, m 3/s, υ is volume of feed, m 3, t is time, s, and C is solute concentration, mol/m 3. C|Z = L is concentration of solute at the outlet of contactor which is inlet for feed tank. Fig. 11 shows the comparisons between simulation results and experimental data reported by Juang and Huang [5]. The module and operational parameters used in the simulations are the same as experimental data. As it can be seen from Fig. 11, simulation results match the experimental data quite well. It confirms that the proposed model is reliable and accurate.
Fig. 11. Comparisons between experimental and simulated results for copper concentrations.
428
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
5. Conclusions A 2D mathematical model was developed for extraction of heavy metals. The solute considered in this work is copper (II). The model is based on solving conservation equations in the membrane contactor. Computational fluid dynamics (CFD) techniques were applied for numerical solution of the equations. Conservation equations of solute including continuity as well as Navier–Stokes equations were derived and solved numerically based on finite element method (FEM). Velocity and concentration distributions of solute were obtained to investigate the effect of process parameters on the extraction of Cu2+. The simulation results showed that decreasing feed velocity increases the solute removal in the contactor. Furthermore, simulation results revealed that membrane contactors are efficient in the extraction of heavy metals. Acknowledgment The authors are grateful of the research council of Islamic Azad University-Farahan Branch for the financial support of this work. References [1] M. Vajda, I. Havalda, R. Macek, Membrane-based solvent extraction and stripping of zinc in a hollow-fibre contactor operating in a circulating mode, Desalination 163 (2004) 19.
[2] L. Sciubba, D.D. Gioia, F. Fava, C. Gostoli, Membrane-based solvent extraction of vanillin in hollow fiber contactors, Desalination 241 (2009) 357. [3] C.B. Patil, P.K. Mohapatra, V.K. Manchanda, Transport of thorium from nitric acid solution by non-dispersive solvent extraction using a hollow fibre contactor, Desalination 232 (2008) 272. [4] B.E. Johnston, Commercial applications of phosphorus-based solvent extraction, Chem. Ind. (London) 20 (1988) 656. [5] R.S. Juang, H.L. Huang, Mechanistic analysis of solvent extraction of heavy metals in membrane contactors, J. Membr. Sci. 213 (2003) 125. [6] R.S. Juang, J.D. Chen, H.C. Huan, Dispersion free membrane extraction: case studies of metal ion and organic acid extraction, J. Membr. Sci. 165 (2000) 59. [7] T.M. Trtic-Petrovic, G.T. Vladisavljevic, M. Tesic, K. Kumric, J.J. Comor, Analysis of concentration boundary layer in thallium (III) extraction with butyl acetate using membrane modules of different length, Desalination 148 (2002) 241. [8] J.J. Guo, C.D. Ho, Theoretical study on membrane extraction of Cu2+ with D2EHPA in laminar flow circular tube modules, Desalination 233 (2008) 247. [9] S. Shirazian, A. Moghadassi, S. Moradi, Numerical simulation of mass transfer in gas–liquid hollow fiber membrane contactors for laminar flow conditions, Simul. Modell. Pract. Theory 17 (2009) 708. [10] A. Gabelman, S.T. Hwang, Hollow fiber membrane contactors, J. Membr. Sci. 159 (1999) 61. [11] R. Prasad, K.K. Sirkar, Hollow fibers solvent extraction: performances and design, J. Membr. Sci. 50 (1999) 153. [12] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, JohnWiley & Sons, 1960. [13] Chemical Engineering Module, COMSOL manual, COMSOL Multiphysics version 3.2. [14] W. Riedl, T. Raiser, Membrane-supported extraction of biomolecules with aqueous two-phase systems, Desalination 224 (2008) 160. [15] M.J. Gonzalez-Munoz, S. Luque, J.R. Alvarez, J. Coca, Simulation of integrated extraction and stripping processes using membrane contactors, Desalination 163 (2004) 1.
Contents lists available at SciVerse ScienceDirect
Desalination journal homepage: www.elsevier.com/locate/desal
Simulation of heavy metal extraction in membrane contactors using computational fluid dynamics Azam Marjani ⁎, Saeed Shirazian Islamic Azad University, Farahan Branch, Department of Chemical Engineering, Farahan, Iran
a r t i c l e
i n f o
Article history: Received 19 March 2011 Received in revised form 12 August 2011 Accepted 13 August 2011 Available online 9 September 2011 Keywords: Heavy metals Extraction Numerical simulation Computational fluid dynamics
a b s t r a c t A new approach based on momentum and mass transfer for modeling of wastewater treatment was developed in this study. The model is based on solving governing equations for solutes using finite element method (FEM). Wastewater treatment medium considered in this study is a microporous membrane contactor made of polypropylene which is a hydrophobic polymer. Modeling was carried out for separation of water contaminants. With solving conservation equations, velocity and concentration distributions were obtained for the solute in the membrane contactor. Modeling findings showed high removal efficiency for copper extraction. The effect of process parameters on the extraction of copper was investigated. Furthermore, simulation results confirmed that membrane contactors are efficient in the removal of water contaminants. © 2011 Elsevier B.V. All rights reserved.
1. Introduction Nowadays, development of clean technologies for removal of heavy metals from water and wastewaters has been of great importance. From environmental perspectives, it is of vital importance to remove heavy metals from wastewaters to prevent environmental problems. Recently, a non-dispersive extraction process utilizing a porous membrane as contactor has attracted many attentions. In this type of process, porous membrane is used as a contact medium between aqueous and organic phases. Membrane contactor is constituted of two parts, i.e. lumen and shell side. Both aqueous and organic solutions flow continuously, the first one in the lumen of fibers and the second in the shell side. Both phases make interface through pores of the fiber wall [1]. Membrane contactors as extraction devices have become a subject of great interest. In these contactors, the membrane acts as a physical barrier between the two phases, i.e. feed and solvent. Membrane contactors provide high surface/volume ratio. In addition, in membrane contactors, the interfacial mass transfer area is relatively high, whereas the hydrophilic or hydrophobic nature of the membrane determines the position of the interface between the feed and the solvent. A major part of the interest towards membrane contactors is due to their capability in setting a dispersion free contact. In addition, the velocities of both phases can be chosen independently, while neither flooding nor unloading problems arise [2,3].
⁎ Corresponding author. E-mail address: [email protected] (A. Marjani). 0011-9164/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.desal.2011.08.032
Extraction of heavy metals using organophosphorus extractants is an important process in hydrometallurgical industries [4]. Conventional solvent extraction processes are carried out using packed towers, mixer-settlers, etc. which try to maximize the contact area of two phases for separation operation. These contactors have some disadvantages such as formation of stable emulsions, thereby inhibiting the phase separation and product recovery; avoid using liquids having similar densities and foaming. Additional limitations present in packed towers include loading requirements and flooding restrictions [5]. Membrane contactors can overcome the drawbacks of conventional solvent extraction processes. Therefore, extraction of heavy metals can be accomplished using membrane contactors. Some researchers have examined the applicability of these contactors for Meal extraction. Juang et al. [5,6] investigated the extraction of copper (II) using di(2-ethylhexyl)phosphoric acid (D2EHPA) in a flat sheet and hollow-fiber membrane contactor. Their results proved the applicability of membrane contactors for removal of copper (II). They also developed a mass transfer model for prediction of membrane extractor performance. Their model was based on resistancein-series model. Trtic-Petrovic et al. [7] used membrane-based solvent extraction for the separation of thallium (III) from chloride-containing acidic solutions with butyl-acetate. They developed a mass transfer model and showed that the existence of concentration boundary layer (CBL) in the aqueous phase inside the hollow fibers is the main limiting factor in TI(III) transport through the interfacial area. Recently Guo and Ho [8] developed a mass transfer model for simulation of Cu 2+ extraction in a membrane contactor. They solved continuity equation for transport of Cu 2+ across the membrane. They
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
used analytical method to solve the governing equation. Parabolic velocity profile for velocity in the contactor was assumed and axial diffusion was neglected which is not a good approximation especially for low flow rates. There are two approaches for modeling of metals extraction in membrane contactors. The first approach is based on resistances-inseries model. This model considers three mass transfer resistances in series for transport of metal from feed to solvent phase. The mass transfer coefficients for all phases are obtained experimentally using correlations which are not very accurate. The second approach is based on solving the conservation equations including mass, momentum and energy for metal in the three sections of membrane contactor, i.e. lumen, membrane and shell side. For this approach numerical method is applied for solving the equations to predict the performance of membrane extractor. There is a definite need for a comprehensive model that can provide a general simulation of the extraction of heavy metals in membrane contactor [9–11]. The main purpose of the current study is to develop and solve a mathematical model for prediction of metal extraction in membrane contactors. The governing equations are then solved numerically using CFD techniques. The velocity field is determined using Navier– Stokes equations. The aim of this simulation is to predict the performance of membrane contactor for extraction of copper.
2. Theory 2.1. Model equations Fig. 1 shows a tubular membrane contactor for extraction of copper. The feed phase containing aqueous solution of Cu 2+ flows inside the inner tube (lumen) while the extractant solution flows countercurrently inside the shell side. Non-wetted mode is assumed for numerical simulation. It is assumed that the extractant-containing phase penetrates the membrane pores and wets them due to
423
hydrophobicity of the membrane. An interface is established at the pores' mouth adjacent to the lumen side. At the interface, chemical reaction occurs between metal ion (Cu 2+) and organometallic complex is formed and diffuses through the membrane pores and is transferred to the shell side. Therefore, the driving force in this process is concentration difference of solute. The mass transfer model is built considering the following assumptions: (1) (2) (3) (4) (5)
Steady state and isothermal conditions Fully developed parabolic velocity profile inside fiber Laminar flow for two phases in the contactor The membrane is completely extraction-phase filled The chemical reaction is very fast and achieves equilibrium.
The main equation that describes the transfer of solute (Cu 2+) from aqueous phase to organic phase is continuity equation. The differential form of continuity equation for solute may be written as [12]: ∂Ci þ ∇⋅ð−Di ∇Ci þ Ci V Þ ¼ Ri ∂t
ð1Þ
where Ci denotes the concentration of solute (mol/m 3), Di denotes its diffusion coefficient (m 2/s), V the velocity vector (m/s) and Ri denotes the reaction term (mol/m 3.s). The velocity vector can be expressed analytically or obtained by coupling a momentum balance to the equation system. This equation is also called convection and diffusion equation. The diffusion coefficient for the dissolved species accounts exclusively for the interaction between the solute and the solvent. In the Maxwell–Stefan convection–diffusion application mode, the interaction between the different dissolved species is also taken into account. Eq. (1) is the main equation of mass transfer. This equation should be solved numerically to obtain the concentration distribution of solute in the membrane contactor. To solve the continuity equation Eq. (1), velocity distribution is needed. Velocity distribution is obtained by solving the momentum equation, i.e. Navier–Stokes equations. Therefore, the momentum and continuity equation should be coupled and solved simultaneously to obtain the concentration distribution of solute. The Navier–Stokes equations describe flow in viscous fluids through momentum balances for each of the components [12]: ∂V T −∇·η ∇V þ ð∇V Þ þ ρðV:∇ÞV þ ∇p ¼ F ∂t ∇:V ¼ 0
ρ
ð2Þ
where η denotes the dynamic viscosity (kg/m.s), V the velocity vector (m/s), ρ the density of the fluid (kg/m 3), p the pressure (Pa) and F is a body force term (N). The continuity equation Eq. (1) is developed for all subdomains of the membrane contactor, i.e. lumen, membrane and shell sides. It may be written as: " # ∂2 Ci 1 ∂Ci ∂2 Ci ∂C Di þ þ 2 ¼ Vz i−shell : r ∂r ∂z ∂r 2 ∂z
ð3Þ
The velocity distribution in the lumen side is assumed to follow Newtonian laminar flow [12]: 2 r VZ−tube ¼ 2u 1− r1
Fig. 1. Tubular membrane contactor for extraction of Cu2+.
ð4Þ
where u is average velocity in the lumen side and r1 is inner radius of lumen (see Fig. 2).
424
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
2.2. Numerical solution of model equations
Fig. 2. Convergence behavior for numerical simulation.
The velocity in the shell side is characterized using Navier–Stokes equations. Eq. (3) should be scaled because of large difference between membrane radius and length. The equation is scaled in the z-direction. A scaling factor is defined as: δ ¼ z=S
ð5Þ
where S is scale factor. Applying the scale factor results in: " # ∂2 Ci 1 ∂Ci ∂2 Ci ∂C Di þ þ ¼ Vz i−shell : r ∂r S2 ∂δ2 S∂δ ∂r 2
ð6Þ
Arranging Eq. (6) results in:
SDi
∂2 Ci 1 ∂Ci Di ∂2 Ci ∂C þ SDi ¼ Vz i−shell : þ 2 r ∂r S ∂δ2 ∂δ ∂r
ð7Þ
Therefore, anisotropic diffusivity is needed for numerical simulation. The diffusivity in r-direction is SDi while the diffusivity in z-direction is Di S. Boundary conditions for numerical simulation are listed in Table 1. Where m is partition coefficient of solute between aqueous phase and organic phase. C1, C2, and C3 are concentration of metal in the lumen, membrane, and shell side respectively. The boundary condition for the convective flux assumes that the mass passing through this boundary is convection-dominated. In other words, it assumes that the mass flux due to diffusion across this boundary is zero. Diffusion mechanism (V = 0) is assumed for mass transfer of solute within the membrane. The chemical reaction for all sections is assumed to be zero (Ri = 0).
Table 1 Boundary conditions of mass transfer equations. Position
Lumen
Membrane
Shell
z = 0 (Inlet) z = L (Outlet) r = r1 r = r2 r = r3
Ci = C0 Convective flux C1 = C2/m – –
Insulated Insulated C2 = C1*m C2 = C3 –
Convective flux Ci = 0 – C3 = C2 ∂Ci ¼ 0 (symmetry) ∂r
The model equations with the appropriate boundary conditions were solved numerically using COMSOL software. This software employs finite element method (FEM) for numerical solutions of model equations. The finite element analysis is combined with adaptive meshing and error control using numerical solver of UMFPACK [9]. This solver is well suited for solving stiff and non-stiff non-linear boundary value problems [9]. An IBM-PC-Pentium4 (CPU speed is 2800 MHz) was used to solve the set of equations. The computational time for solving the set of equations was about 2 min. Table 2 shows parameters used in the numerical simulation. The maximum number of iterations was set at 25 and the relative tolerance was 1E−6. After about 9 iterations the system reaches convergence. Fig. 2 shows the convergence behavior for numerical simulation. Fig. 2 implies that the system is not stiff and provides confidence in the solution process and ensures consistency. Fig. 3 shows a segment of the mesh used to determine the solute transport behavior in the membrane contactor (MC). It should be pointed out that the COMSOL mesh generator creates triangular meshes that are isotropic in size. A large number of elements are then created with scaling. A scaling factor of 120 has been employed in z-direction due to large difference between r and z. COMSOL automatically scales back the geometry after meshing. This generates an anisotropic mesh around 4397 elements. Adaptive mesh refinement in COMSOL, which generates the best and minimal meshes, was used to mesh the contactor geometry. For numerical simulation, it is important to select a mesh that minimizes the error in the calculations. It is not easy to specify mesh sizes manually to minimize the error. The algorithm must resolve the solution in great detail only on small portions of the domain. Adaptive mesh generation identifies the regions where high resolution is needed and produces an appropriate mesh. The adaptive solver works on linear or nonlinear stationary PDE problems (including parameter-dependent problems) or eigenvalue PDE problems using the coefficient or general solution forms. The adaptive mesh refinement performs the following iterative algorithm for numerical simulation of model's equations: 1. Solve the problem on the existing mesh using the stationary solver. 2. Evaluate the residual of the PDE on all mesh elements. 3. Estimate the error in the solution on all mesh elements. The computed error estimate is really just an error indicator because the estimate involves an unknown constant. 4. Terminate execution if it has made the requested number of refinements or if it has exceeded the maximum number of elements. 5. Refine a subset of the elements based on the sizes of the local error indicators. 6. Repeat from Step 1 [13].
Table 2 Parameters used in the simulations [5,6]. Extractor length (cm) Fiber o.d. (mm) Fiber i.d. (mm) Module i.d. (mm) Scale factor Kex Porosity (ε) Tortuosity (τ) DCu2+-water (m2/s) Dcomplex-organic (m2/s) Dcomplex-membrane (m2/s) Extractant concentration (mol/m3) Inlet concentration of solute (mol/m3)
30 20 22 30 120 4.4E−4 0.67 1.5 3.25E−8 2E−10 Dcomplex-organic (ε/τ) 30 1
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
425
(extractant) flows from the other side (z = L) where the concentration of solute is assumed to be zero. As the feed flows through the lumen side, solute (Cu 2+) is transferred towards the membrane due to the concentration difference. At the membrane-aqueous phase interface chemical reaction occurs and metallic complex is formed. The formed complex diffuses through the membrane pores which is filled with the organic phase and reaches at the shell side. Finally, the complex is swept by the moving extractant and leaves the extractor. Fig. 4 confirms that the concentration decreases along the extractor because of the extraction of solute. Fig. 4 also indicates that increasing feed flow rate decreases the solute removal in the contactor. Increasing feed flow rate decreases residence time of feed phase in the contactor which in turn results in lower mass transfer of solute [14,15]. The extraction efficiency of solute which is defined the ratio of the solute transfer from the feed phase to extraction phase to total solute in the initial feed phase is calculated using numerical simulation. The extraction efficiency for solute is calculated 85.8% for extraction of Cu 2+ using membrane extractor. 3.2. Contours of concentration in the lumen side
Fig. 3. Magnified segment of the mesh used in the numerical simulation (There are 4397 elements in total for the whole MC domain; z-direction scale factor= 120; the three domains from left to right are lumen side, membrane and shell side, respectively).
To investigate the concentration of solute in the lumen side of the membrane contactor, contours of solute concentration is calculated. Fig. 5 shows the contours of solute concentration in the lumen side where the feed phase flows with laminar flow. As it can be seen from Fig. 5, concentration gradient is significant at the regions near
3. Results and discussion 3.1. Concentration profile of solute along the membrane contactor Dimensionless concentration profile (C/C0) of solute is shown in Fig. 4. The concentration profile is obtained in the lumen side of the membrane contactor where the aqueous feed flows. It should be pointed out that the flow pattern in the contactor (extractor) is parallel and counter current. The feed phase including aqueous solution of Cu 2+ flows from one side of the contactor (z = 0) where the concentration of solute is the highest (C0), whereas the organic phase
Fig. 4. Axial concentration distribution of solute along the extractor. Organic flow rate = 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature = 298 K, Pressure = 1 atm.
Fig. 5. Contours of concentration in the lumen side. Feed flow rate = 1E−8 m3/s, organic flow rate = 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature= 298 K, Pressure= 1 atm.
426
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
the inlet of lumen side. High mass transfer flux occurs in these regions of lumen because of higher driving force. As the feed flows along the lumen side, it is transferred to the shell side and causes lower concentration gradient at the regions near the lumen outlet. 3.3. Profile of concentration in the shell side Fig. 6 indicates concentration profile of solute in the shell side of contactor. Obviously, concentration difference at a narrow region near the fiber wall is sharp and significant. This is an evidence for the formation of concentration boundary layer in the shell side. At the regions far from the fiber wall, solute concentration is zero. It reveals that the solute cannot be transferred into the distances far from fiber wall and is swept by the solvent flow. Another reason for such behavior could be attributed to this fact that the solubilizing capacity of the extractant is high and partition coefficient of solute between feed and organic phase is relatively high. 3.4. Streamline of solute concentration in the lumen side Streamline of solute concentration in the lumen side of membrane contactor is illustrated in Fig. 7. It is seen that the solute is transferred from the bulk of feed phase towards the membrane in the r-direction because of concentration difference. The mass transfer mechanism in the r-direction is governed by molecular diffusion. It means that the diffusion is the predominant in the radial direction while in axial direction the convection is predominant. It is clearly seen that at the regions near the lumen center, convection in the z-direction is significant and causes mass transfer of solute towards the lumen outlet. However, to increase the removal rate of solute, concentration gradient in r-direction should be enhanced. This causes the mass transfer of solute to membrane and shell side. This can be accomplished by choosing strong and effective extractant. 3.5. Velocity field in the extractor The contours of velocity in the shell side of membrane contactor are illustrated in Fig. 8. Fig. 9 also indicates the profile of organic phase in the shell side at different axial positions. The velocity field in the shell side of membrane contactor was simulated by solving the Navier–Stokes equations. As it can be seen from Figs. 8 and 9, the velocity profile is almost parabolic with a
Fig. 7. Streamline of solute concentration in the lumen side. Feed flow rate = 1E −8 m3/s, organic flow rate = 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature = 298 K, Pressure = 1 atm.
mean velocity which increases with the membrane length because of continuous fluid permeation. Fig. 9 also reveals that at the inlet regions in the liquid side, the velocity is not developed. After a short distance from the inlet, the velocity profile is fully developed (see Fig. 9). As observed, the model considers the entry effects on the hydrodynamics of fluid flow in the shell side. Furthermore, Fig. 8 reveals the effect of viscous forces at the region adjacent the shell wall. The effect of viscous force in these regions is significant and causes zero velocity (no slip) at the shell wall. 3.6. Effect of porosity-to-tortuosity ratio In membrane contactors, porosity and tortuosity of membrane is of great importance. These parameters have effect on the extraction efficiency of solute [10]. To investigate the effect of porosity and tortuosity on the performance of process, porosity-to-tortuosity ratio varied between 0.05 and 0.80 has been considered in the simulations. Fig. 10 shows the effect of the porosity-to-tortuosity ratio on the extraction efficiency of copper. The removal of copper has grown about 7 times when the latter ratio increases from 0.05 to 0.80. The effective diffusion coefficient (Dcomplex-membrane) of the membrane is calculated through the porosity and tortuosity of the membrane [10]: Dcomplex−membrane ¼ Dcomplex−organic ðε=τÞ:
Fig. 6. Concentration profile of solute in the shell side. Feed flow rate= 1E−8 m3/s, organic flow rate= 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature= 298 K, Pressure = 1 atm.
ð8Þ
As it can be seen from Eq. (8), the effective diffusion coefficient is a function of porosity and tortuosity of membrane. As the porosity/tortuosity ratio increases, the membrane diffusivity and thus the mass transfer of complex through the membrane increases; i.e. the extraction efficiency of solute increases. Indeed, when the porosity/tortuosity ratio increases, the membrane mass transfer resistance decreases. Therefore, the total resistance to the mass transfer of solute becomes lower.
Extraction efficiency of solute (%)
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
427
90 80 70 60 50 40 30 20 10 0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Porosity / tortuosity ratio Fig. 10. Effect of porosity-to-tortuosity ratio on extraction efficiency of solute. Feed flow rate = 1E−8 m3/s, organic flow rate = 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, extractant concentration = 30 mol/m3, Temperature = 298 K, Pressure = 1 atm.
4. Model validation To validate the modeling findings obtained from mass transfer model developed here, the simulation results were compared with the experimental data obtained from literature for extraction of copper using a hollow-fiber membrane contactor. The experimental data are for recycling mode of operation. An additional equation for feed tank is required for simulation of recycling mode. The mass balance equation for feed tank containing aqueous solution of copper (solute) can be written as follows: υ
Fig. 8. Contours of velocity in the shell side. Feed flow rate = 1E−8 m3/s, organic flow rate = 1E−7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature = 298 K, Pressure = 1 atm.
Fig. 9. Velocity profile in the shell. Feed flow rate = 1E−8 m3/s, organic flow rate = 1E −7 m3/s, inlet concentration of solute = 1 mol/m3, Temperature = 298 K, Pressure = 1 atm.
dC tank ¼ QC jZ¼L −QC tank dt
@t ¼ 0;
Ctank ¼ Cint
ð9Þ ð10Þ
where Qis volumetric flow rate, m 3/s, υ is volume of feed, m 3, t is time, s, and C is solute concentration, mol/m 3. C|Z = L is concentration of solute at the outlet of contactor which is inlet for feed tank. Fig. 11 shows the comparisons between simulation results and experimental data reported by Juang and Huang [5]. The module and operational parameters used in the simulations are the same as experimental data. As it can be seen from Fig. 11, simulation results match the experimental data quite well. It confirms that the proposed model is reliable and accurate.
Fig. 11. Comparisons between experimental and simulated results for copper concentrations.
428
A. Marjani, S. Shirazian / Desalination 281 (2011) 422–428
5. Conclusions A 2D mathematical model was developed for extraction of heavy metals. The solute considered in this work is copper (II). The model is based on solving conservation equations in the membrane contactor. Computational fluid dynamics (CFD) techniques were applied for numerical solution of the equations. Conservation equations of solute including continuity as well as Navier–Stokes equations were derived and solved numerically based on finite element method (FEM). Velocity and concentration distributions of solute were obtained to investigate the effect of process parameters on the extraction of Cu2+. The simulation results showed that decreasing feed velocity increases the solute removal in the contactor. Furthermore, simulation results revealed that membrane contactors are efficient in the extraction of heavy metals. Acknowledgment The authors are grateful of the research council of Islamic Azad University-Farahan Branch for the financial support of this work. References [1] M. Vajda, I. Havalda, R. Macek, Membrane-based solvent extraction and stripping of zinc in a hollow-fibre contactor operating in a circulating mode, Desalination 163 (2004) 19.
[2] L. Sciubba, D.D. Gioia, F. Fava, C. Gostoli, Membrane-based solvent extraction of vanillin in hollow fiber contactors, Desalination 241 (2009) 357. [3] C.B. Patil, P.K. Mohapatra, V.K. Manchanda, Transport of thorium from nitric acid solution by non-dispersive solvent extraction using a hollow fibre contactor, Desalination 232 (2008) 272. [4] B.E. Johnston, Commercial applications of phosphorus-based solvent extraction, Chem. Ind. (London) 20 (1988) 656. [5] R.S. Juang, H.L. Huang, Mechanistic analysis of solvent extraction of heavy metals in membrane contactors, J. Membr. Sci. 213 (2003) 125. [6] R.S. Juang, J.D. Chen, H.C. Huan, Dispersion free membrane extraction: case studies of metal ion and organic acid extraction, J. Membr. Sci. 165 (2000) 59. [7] T.M. Trtic-Petrovic, G.T. Vladisavljevic, M. Tesic, K. Kumric, J.J. Comor, Analysis of concentration boundary layer in thallium (III) extraction with butyl acetate using membrane modules of different length, Desalination 148 (2002) 241. [8] J.J. Guo, C.D. Ho, Theoretical study on membrane extraction of Cu2+ with D2EHPA in laminar flow circular tube modules, Desalination 233 (2008) 247. [9] S. Shirazian, A. Moghadassi, S. Moradi, Numerical simulation of mass transfer in gas–liquid hollow fiber membrane contactors for laminar flow conditions, Simul. Modell. Pract. Theory 17 (2009) 708. [10] A. Gabelman, S.T. Hwang, Hollow fiber membrane contactors, J. Membr. Sci. 159 (1999) 61. [11] R. Prasad, K.K. Sirkar, Hollow fibers solvent extraction: performances and design, J. Membr. Sci. 50 (1999) 153. [12] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, JohnWiley & Sons, 1960. [13] Chemical Engineering Module, COMSOL manual, COMSOL Multiphysics version 3.2. [14] W. Riedl, T. Raiser, Membrane-supported extraction of biomolecules with aqueous two-phase systems, Desalination 224 (2008) 160. [15] M.J. Gonzalez-Munoz, S. Luque, J.R. Alvarez, J. Coca, Simulation of integrated extraction and stripping processes using membrane contactors, Desalination 163 (2004) 1.