Modeling of mass transport of aqueous solutions of multi-solute organics through reverse osmosis membranes in case of solute–membrane affinity

Journal of Membrane Science 267 (2005) 27–40

Modeling of mass transport of aqueous solutions of multi-solute organics through reverse osmosis membranes in case of solute–membrane affinity Part 1. Model development and simulation夽 H. Mehdizadeh a,∗ , Kh. Molaiee-Nejad a , Y.C. Chong b a b

Chemical Engineering Department, Tarbiat Modarres University, P.O. Box 14115-143, Tehran, Iran School of Chemical Engineering, University of Science of Malaysia, 31750 Tronoh, Perak, Malaysia Received in revised form 7 February 2005; accepted 10 March 2005 Available online 22 July 2005

Abstract A two-dimensional mathematical model is developed for transport of multi-solute liquid (aqueous) solutions through reverse osmosis membranes in which strong affinity may exist between the membrane polymer and the solutes. The final objective of this research study is to predict the performance of reverse osmosis (RO) or nanofiltration (NF) membranes in the case of multi-solute aqueous systems in the presence of strong solute–membrane attraction. Modeling of membrane transport in this case is complex because of the interactions between solute, solvent and membrane. The model, which assumes a micro-porous structure for the membrane, is an extension of a single-solute model based on the preferential sorption—capillary flow mechanism, and takes detailed solute–solvent–membrane interactions into account. The model is believed to be the first to have been able to describe the anomalous behavior of such systems for multi-component aqueous solutions. The developed model is used to simulate the performance of cellulose acetate membrane in aqueous toluene–benzene systems by varying the operating pressure and mole fraction of toluene/benzene system. Generally, the separation of both solutes will decrease as the operating pressure is increased, with toluene having higher separations than benzene. The permeation flux also decreases with increasing operating pressure. However, the permeation flux of mixture increases as solutes concentration in feed solution is increased. The solutes separation values are higher at higher solutes concentrations of the feed. © 2005 Elsevier B.V. All rights reserved. Keywords: Membrane; Reverse osmosis; Modeling; Cellulose acetate; Multi-component; Affinity

1. Introduction The phenomenon of strong solute–membrane affinity has been studied and reported by some researchers [1–16] when certain low-molecular weight organic solutes, such as toluene, benzene, cumene, and phenol derivatives, exist in water. These organic compounds are common contaminants in wastewater which could be generated by a variety of chem夽 Presented at the 16th International Congress of Chemical and Process

Engineering (CHISA 2004), Prague, Czech Republic, 22–26 August 2004, Paper # P3.57, p. 665. ∗ Corresponding author at: Chemical Engineering Department, University of Qatar, P.O. Box 2713, Doha, Qatar. Fax: +974 417 4852101. E-mail address: [email protected] (H. Mehdizadeh). 0376-7388/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.memsci.2005.03.059

ical and petrochemical industries. For example, chlorophenols, known as a dangerous, toxic, environmentally persistent material, are widely used in manufacturing of pesticides, insecticides, and so on. Most of such materials are dangerous toward living organisms and some of them are carcinogenic. The reverse osmosis performance in such systems is markedly different than the normal systems in which solute molecules are rejected by the membrane. Some of the main characteristics behavior are summarized below: the separation of these organic solutes decreases as the operating pressure increases; permeate flux is usually much less than pure water (solvent) flux even when osmotic effects are small; partition coefficients (ratio of concentration in the membrane over the concentration in the neighboring bulk solution at equilibrium) are larger than unity; and separation can be

28

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

positive, zero, or negative, depending on selective operating conditions [1–16] as is discussed in the following sections. Modeling of reverse osmosis can be carried out by three approaches: (i) based on black box irreversible thermodynamics, IT, such as the “Spiegler–Kedem model”, (ii) based on nonporous membranes, such as the “solution–diffusion” (SD) model, and (iii) based on porous membranes, such as the “Kimura–Sourirajan analysis”, KSA, and the “finely porous model”, FPM [17]. Some of these models have identical mathematical equations despite their contradictory premises (such as the SD model and the KSA, or the FPM relationship and the IT Kedem–Spiegler model). In the existing theories, the two main viewpoints are models based on the solution–diffusion mechanism as proposed by Lonsdale et al. [18] and models based on preferential sorption-capillary flow (PSCF) mechanism as proposed by Sourirajan et al. [17,19]. Despite the shortcomings noticed for each transport model, both sides have been able to predict the membrane systems to some extent of success. For example, the SD model and the KSA both predict a 100% separation at very large permeation fluxes, which contradicts some experimental data for some aqueous systems with inorganic solutes; however these models are still used in many cases due to their simplicity. A two-dimensional transport model, based on the PSCF mechanism which takes into account the strong affinity which may exist between solute molecules and membrane polymer in a single-solute liquid system has been presented by Mehdizadeh and Dickson [20]. In the present work, the single-solute model is extended to the more practical case of multi-solute systems as discussed in the following Theory section. The single-solute model [20] is based on an original model developed for the normal RO case of solute repulsion by membrane [21]. The original model [21] has been evaluated for different normal systems which has proved to be efficient and promising in prediction of membrane performance (see for example [22–23]).

2. Theory Modeling and simulation of membrane separation systems are done for a variety of reasons the most important of all to get a reliable prediction of the membrane system performance and a reliable base for design of membrane modules. In simulation studies, model parameters are assumed to be known, and the behavior of an existing membrane system (at least conceptually) is determined under certain operating conditions by applying the model with its known parameters. On the other hand, if the model parameters are not known (or cannot be determined a priori) experimental data can be used to obtain the best values of the parameters by using the transport model in conjunction with an optimization code. In either case (simulation or optimization), the preferred model, among a number of existing models, is the one whose parameters are less and not dependent on the driving forces involved in the membrane process, with

the latter characteristic as the most important feature of the model. For example, the membrane models based on irreversible thermodynamics are highly dependent on the driving forces in reverse osmosis (which are pressure and concentration gradients), probably due to the application of Onsager reciprocal relations in deriving the model equations; and this disadvantage highly restricts the practical applications of the model. The phenomenon of solute–membrane affinity has been noticed by some researchers to happen when certain lowmolecular weight organic solutes exist in water, such as toluene, benzene, cumene, and phenol derivatives. The reverse osmosis performance in such systems is markedly different than the usual systems in which solute is rejected by the membrane. The anomalous behavior has been attributed to the non-polar or polar interactions between the organic molecules and the backbone of the polymer material. The characteristic behavior for such systems can be summarized as: increasing the operating pressure decreases the separation (unlike the usual case in solvent–membrane affinity systems); permeate flux is usually much less than the pure water (solvent) flux even when osmotic effects are small; partition coefficients are larger than unity; and separation can be positive, zero, or negative (see for example [4–12,16]). Lonsdale et al. [1] first reported such solute–membrane affinity behavior for phenol–water separation with cellulose acetate membranes, but they could not model the data. Since then, several researchers modeled the data using the fourparameter FPM or IT relationship (e.g. [24]). However, none of the models, except the single-solute model by Mehdizadeh and Dickson [20], has been able to predict the behavior that separation can decrease from positive to negative values with increasing applied pressure. However, the successful model [20] has been derived for single-solute systems, and the real world deals with multi-solute systems, such as the wastewater treatment plants. In the present article, the single-solute model for the affinity system [20] is extended to the more practical case of multi-solute systems. The derivation of the governing equations is given in the following section. Finally, simulation studies are presented based on the developed model. Experimental evaluation of the new model will be presented in a second paper. 2.1. Mathematical modeling In this section, the mathematical model for a multi-solute system, with strong solute–membrane affinity, is derived in a manner paralleling the single-solute model [20]. Starting with the radial component of molar flux for each solute, and then the axial component of the flux together with the equation of continuity, equations of concentration profile and permeation flux, for each solute, are derived. Finally, following the derivation of the differential equation for the fluid velocity inside the membrane, the equations required to describe the overall flux and separation are presented.

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

29

Eqs. (2) and (3) can be used to define the partition coefficients on the high- and low-pressure sides of the membrane:
 Ci (ρ, ξ = 0) ϕi (r, 0) = exp − K2i (ρ) = Ci2 RT = exp[−Φi (ρ, 0)]

(4)


 Ci (ρ, ξ = 1) ϕi (r, τ) K3i (ρ) = = exp − Ci3 RT = exp[−Φi (ρ, 1)]

(5)

where the dimensionless potential function: Fig. 1. Cylindrical coordinate system in a membrane pore; the left end of the pore faces feed solution, and the right end faces permeate solution [20].

Consider a liquid feed solution consisting of n solutes plus a solvent (e.g., water), in which the solutes show strong affinity towards a microporous membrane when they come into close contact with the membrane polymer. The membrane is modeled as microporous with pores modeled as perfect cylinders of radius RW and length τ, as shown in Fig. 1. Only the skin layer of the membrane is shown in Fig. 1 with a length equal to τ which takes tortuosity into account, and tortuosity is defined as the ratio of the actual distance a molecule travels between two ends of a membrane to the shortest distance (thickness) between the two ends. 2.1.1. Radial component of solute flux Since τ  RW , it is fairly reasonable to assume that the radial component of solute flux is almost zero. The assumption of radial equilibrium has been discussed by a few researchers (e.g. [25]). Therefore, radial equilibrium states that for solute i in the pore: Ji,r ∂Ci (r, z) Ci (r, z) ∂ϕi (r, z) = + = 0, −DiB ∂r RT ∂r

i = 1, n

(1)

where the pore wall potential function on solute i, ϕi (r, z), which is the source for the net body force acting on solute i by the pore wall, is negative when the solute i is attracted by the pore wall (and it is positive if the solute is repelled by the pore wall). Integration of Eq. (1) leads to the following Boltzmann equations at the pore entrance (z = 0):
 ϕi (r, 0) Ci (r, 0) = Ci2 exp − , i = 1, n (2) RT and the pore exit (z = τ):
 ϕi (r, τ) , Ci (r, τ) = Ci3 exp − RT

i = 1, n

(3)

where Ci2 and Ci3 are the solute i concentrations at the membrane–feed side and at the permeate side, respectively. Derivations of Eqs. (2) and (3) parallel that in Appendix A of a previous paper [21].

ϕi (r, z) , i = 1, n (6) RT is written in terms of the dimensionless coordinate system ρ and ξ: r ρ= (7) RW z (8) ξ= τ Therefore, according to Eqs. (4) and (5), the partition coefficients are different if the potential functions are different at the ends of a pore. Effective (average) partition coefficients for solute i (K2i and K3i ) can be obtained by averaging Eqs. (4) and (5) over the pore cross-sectional area:  1−λi  1−λi K2i (ρ)ρ dρ 0 K2i = e−Φi (ρ,0) ρ dρ (9) =2 1 0 ρ dρ 0  1−λi  1−λi K3i (ρ)ρ dρ 0 K3i = e−Φi (ρ,1) ρ dρ =2 (10) 1 0 ρ dρ 0 Φi (ρ, ξ) =

where λi is the ratio of molecular radius of solute i to the pore radius. For the imaginary case of zero potential field, Eqs. (9) and (10) reduce to values for size exclusion by steric hindrance: K2i = K3i = (1 − λi )2

(11)

which is the case in simple filtration or size screening. 2.1.2. Axial component of solute flux and equation of continuity Four forces acting on solute i inside a pore are: the total driving force of solute i, Fi (r, z), the frictional force between the solvent B and the solute i, FiB (r, z), the summation of  all frictional forces between other solutes and solute i, nj=1(j=i) Fij (r, z), and the frictional force between the solute i and the pore wall, FiM (r, z). The total driving force of solute i is given by the gradient of the solute chemical potential inside the pore [25], Fi (r, z) = −

∂µi (r, z) , ∂z

i = 1, n

(12)

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

30

where the total chemical potential is given by (see for example [26]): µi = µ0i + RT ln ai + vi P + ϕi (r, z),

i = 1, n

(13)

The vi P term accounts for the transport of solute molecules by pressure forces, and the ϕi term takes into account the effect of pore wall potential on the solute transport. By assuming an ideal solution, the gradient can be written for the isothermal system as ∂µi (r, z) 1 ∂[RTCi (r, z)] ∂P(r, z) = + vi ∂z Ci (r, z) ∂z ∂z +

∂ϕi (r, z) ∂z

(14)

which is responsible for driving the solute i through the pore. Then, using Eq. (12), upon using van’t Hoff’s equation for solution osmotic pressure, the driving force becomes: RT ∂πi (r, z) ∂P(r, z) ∂ϕi (r, z) Fi (r, z) = − − vi − πi (r, z) ∂z ∂z ∂z

(15)

To derive the axial component of the solute flux, Eq. (15) is used in a force balance on the solute inside the pore, as presented in Appendix I. The result is:  πi (ρ, ξ) ∂P(ρ, ξ) 1 DiB ∂πi (ρ, ξ) Ji,z (ρ) = − vi − τRT bi (ρ) ∂ξ RT ∂ξ  ¯ D + α(ρ)πi (ρ, ξ) (16) DiB where bi (ρ) and α(ρ) are the friction function and fluid velocity defined as (both are dimensionless):  χiB + χiM (ρ) + nj=1(=i) χij DiB (17) = bi (ρ) = DiM (ρ) χiB α(ρ) =

uB (ρ)τ ¯ D

(18)

¯ is defined, in this paper, as An average solute diffusivity, D, n j=1 DjB ¯ = D (19)

n−1 n D j=1 jB The equation of continuity of solute i states that: ∇ · Ji = 0

Assuming that pressure changes linearly in ξ [21], then: ∂P(ρ, ξ) = P(ρ, 1) − P(ρ, 0) = −P ∂ξ  n

n + πi2 σi2 (ρ) − πi3 σi3 (ρ) i=1

(23)

i=1

which is a result of the radial equilibrium assumption (see [25]). σ i2 (ρ) and σ i3 (ρ) are the local reflection (Staverman) coefficients, for each solute i, at the two sides of the membrane, defined as σi2 (ρ) = 1 − e−Φi (ρ,0)

(24)

−Φi (ρ,1)

(25)

σi3 (ρ) = 1 − e

Therefore, Eq. (22) is written as ∂πi (ρ, ξ) ∂2 πi (ρ, ξ) − {ωi (ρ) + mi α(ρ)} =0 ∂ξ 2 ∂ξ

(26)

where mi is defined as mi =

¯ D DiB

(27)

and ωi represents the solute velocity induced by the net pressure forces:    n n vi πi2 σi2 (ρ) − πi3 σi3 (ρ) P − ωi (ρ) = RT i=1 i=1 (28) Therefore, larger values of partial molar volume or operating pressure will make larger pressure-induced transport of the solute. Then, the differential equation (26) can be solved subject to the boundary conditions (similar to Eqs. (2) and (3)): πi (ρ, 0) = π2 e−Φi (ρ,0)

(29)

−Φi (ρ,1)

(30)

πi (ρ, 1) = π3 e

to give the profiles of osmotic pressure or concentration, for solute i, inside the pore:   e[mi α(ρ)+ωi (ρ)]ξ − 1 πi (ρ, ξ) = πi2 − [πi2 − πi3 Ki∗ (ρ)] [m α(ρ)+ω (ρ)] i e i −1

(20)

× e−Φi (ρ,0)

(31)

which reduces to: ∂Ji,z =0 ∂ξ

(21)

Then, substituting for the flux, from Eq. (16), into the above equation will yield:   ∂2 πi (ρ, ξ) vi ∂ ∂P(ρ, ξ) π + (ρ, ξ) i ∂ξ 2 RT ∂ξ ∂ξ ¯ ∂πi (ρ, ξ) D − α(ρ) =0 DiB ∂ξ

(22)

 Ci (ρ, ξ) =

Ci2 − [Ci2 − Ci3 Ki∗ (ρ)]

×e−Φi (ρ,0)

e[mi α(ρ)+ωi (ρ)]ξ − 1 e[mi α(ρ)+ωi (ρ)] − 1



(32)

where Ki∗ (ρ) is the ratio of the local partition coefficients at the two ends of the pore: Ki∗ (ρ) =

K3i (ρ) e−Φi (ρ,1) = −Φ (ρ,0) K2i (ρ) e i

(33)

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

This ratio can be small or large which, in turn, leads to positive or negative separations. It is mainly dependent on the potential function values at the two ends of a pore. Finally, substituting from Eq. (31) into Eq. (16), and then combining with Eq. (23), explicit expression for the solute flux can be derived as   DiB mi α(ρ) + ωi (ρ) πi2 − Ki∗ (ρ)πi3 Ji,z (ρ) = πi2 + [m α(ρ)+ω (ρ)] i bi (ρ) τRT e i −1 × e−Φi (ρ,0)

(34)

2.1.3. Derivation of velocity profile As derived in Appendix II, the following nonlinear differential equation for the fluid velocity can be derived by applying a force balance in the z-direction on the fluid element in the annular region between z and z + dz and between r and r + dr (see Fig. 1):   2  n ∂ α(ρ) 1 ∂α(ρ) 1  P 1 + + − [πj2 σj2 (ρ) ∂ρ2 β1  π 2 ρ ∂ρ π2 j=1

 

 n  1 1 − πj3 σj3 (ρ)] − 1− [mj α(ρ) + ωj (ρ)]  β1 bj (ρ)  × 1+

j=1

πj2 − Kj∗ (ρ)πj3 πj2 (e[mj α(ρ)+ωj (ρ)] − 1)



π2 =

n

concentration effects and interactions between solutes, solvent, and membrane, in a nonlinear way. These effects will finally lead to reduced permeation flux with solutes present. 2.1.4. Solvent and solute flux through the membrane In order to determine the fluxes over the whole surface area of a membrane, the solvent and solute fluxes are initially derived for flow through a single pore as functions of radial position, and then the relationships are integrated over the surface area of the single pore to find the average fluxes. Finally, the fluxes for a single pore are related to the permeation fluxes through the membrane by using the fractional pore area, ε. The average velocity of the solvent in a single pore is: 1   1 ¯ uB (ρ)ρ dρ D α(ρ)ρ dρ (40) =2 uB = 0  1 τ 0 ρ dρ 0

and the average solvent flux is: JB = uB C = where  l1 =

e−Φj (ρ,0) = 0

(35)

(36)

(37)

¯ σ i2 (ρ), σ i3 (ρ), ωi (ρ), and Ki∗ (r) are defined and bi (ρ), α(ρ), D, by Eqs. (17)–(19), (24), (25), (28) and (33), respectively. The boundary conditions for the velocity profile, α(ρ), are:

dα(ρ) =0 dρ

at ρ = 1 at ρ = 0

1

α(ρ)ρ dρ

(42)

(38) (39)

In the special case of single-solute system (n = 1), negligible νi , and only radial dependent potential function Φ(ρ) [that is, K2 (ρ) = K3 (ρ)], the velocity profile, Eq. (35), correctly reduces to the form in the original model (Eq. (26) in [21]). Furthermore, Eq. (35) reduces to the Poiseuille parabolic velocity profile in the imaginary case of no potential function (Φ = 0), free solution diffusivity (b = 1), and pure solvent flow (π2 = π3 = 0). The second and third (last) terms in brackets, in Eq. (35), count for distortion of the velocity profile from the Poiseuille profile as a result of strong

2DiB l2i τRT

where  l2i =

1−h¯

0

j=1

α(ρ) = 0

(41)

The average flux through the pore for solute i is given by integrating Eq. (34) over the pore area: Ji =

πj2

¯ 2CD l1 τ

0

where ¯ ηD β1 = 2 RW π2

31

(43)

mi α(ρ) + ωi (ρ) bi (ρ)

 πi2 +

× e−Φi (ρ,0) ρ dρ

πi2 − Ki∗ (ρ)πi3 [m e i α(ρ)+ωi (ρ)] − 1



(44)

and h¯ is the dimensionless thickness of the adsorbed layer of the solutes on the pore wall: h h¯ = RW

(45)

Then, the total solute flux will be: n j=1

1 DjB l2j τRT n

Jj =

(46)

j=1

Using the fractional pore area, ε, the above average fluxes can be related to the fluxes through the membrane, for solvent: NB = εJB

(47a)

which can be written, as shown elsewhere [21], for pure water flux through the membrane as

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

32

CR2W P 8η(τ/ε)

NP =

(47b)

and for total solutes: n

NS =

n Nj = ε J j

j=1

(48)

j=1

Therefore, the total permeation flux through the membrane, NT , is:  ¯ 1 2 nj=1 DjB l2j 2CDl NT = (NS + NB ) = + (49a) RT (τ/ε) τ/ε Dividing Eq. (49a) by Eq. (47b) will yield the flux ratio as   n 16η  1 NT ¯ 1 = 2 (49b) DjB l2j + Dl NP RW P CRT j=1

On the other hand, each solute concentration in the permeate can be determined from a material balance on each solute flux across the membrane,     Ni Ji =C (50) Ci3 = C n NS + N B j=1 Jj + JB or, using Eqs. (44)–(49), Ci3 =

l2i  , (CRT )l1 + nj=1 dj l2j

i = 1, n

(51)

where dj =

DjB ¯ D

(52)

Eq. (51) presents a practical way to compute permeate concentration for each solute if the velocity profile is known from Eq. (35), so that the total permeate concentration can be found using the following equation: C3 =

n

(53)

Cj3

j=1

with Ci3 given by Eq. (51). Finally, the separation achieved by a membrane, for each solute, is computable from: fi =

Ci2 − Ci3 , Ci2

i = 1, n

(54)

The adsorbed layer thickness, h, can also be determined from the relationship between the total solution flux, NT , and pure solvent flux, NP , as given by Thiel et al. [27]:   NT h 4 = 1− (55) NP RW 2.1.5. Potential and friction functions It is desirable to estimate the potential function, Φ, from the basic physicochemical characteristics of a membrane, such as the surface-charge density and electrochemical double-layer forces (for the case of electrolytic solutions) or the Hamaker constant and van der Walls forces (for the case of uncharged organic solutes). However, unless such clear relationships for the potential function are known, in advance, empirical relationships are unavoidable due to the complexity of such relationships. In the original model, the following form of the potential function was used for the case of single-solute repulsion by membrane [21] with positive values (i.e., repulsion) for the potential parameter, θ 1 ,   θ1 eρ2 /2 when ρ < 1 − λ Φ(ρ) = RW (56)  ∞ when ρ ≥ 1 − λ The following form was proposed for the case of single-solute attraction by membrane [20]:   θi eρ2 /2 F {C , C , ξ} when ρ < 1 − λ i i2 i3 i Φi (ρ, ξ) = RW  ∞ when ρ ≥ 1 − λi (57) with negative values (i.e., attraction) for θ i , where Fi is a function whose value depends on the solute concentration inside the pore, Ci (ρ, ξ). However, Ci (ρ, ξ) needs to be determined from Eq. (32), which in turn depends on Φi values. Therefore, the equations are coupled and not easy to solve. In order to avoid such difficulty, and realizing that the equation for solute concentration is already coupled with equation of fluid velocity, the following relationship for Fi is proposed so that the whole set of governing equations can be solved in a manageable way: Fi {Ci2 , Ci3 , ξ} = eγi ξ {1−γi [ξ ln Ci3 + (1 − ξ) ln Ci2 ]} (58) where γ i is the additional potential parameter which adjusts the variation of the concentration profile. The above function has been obtained by trial and error, in the present work, and has shown to be the best representing, yet simple, function so far. Therefore for a multi-solute system, where each type of solute may experience a different potential field, the proposed equation for the potential function is:

  θi e(γi ξ+ρ2 /2) {1 − γ [ξ ln C + (1 − ξ) ln C ]} when ρ < 1 − λ i i i3 i2 Φi (ρ, ξ) = RW  ∞ when ρ ≥ 1 − λi

(59)

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

Two important values, frequently needed in the present model, are the values of the potential functions at the pore entrance and exit; they are explicitly derived from Eq. (59) as Φi (ρ, 0) =

θi ρ2 /2 e {1 − γi ln Ci2 } RW

when ρ < 1 − λi (60)

Φi (ρ, 1) =

θi (γi +ρ2 /2) e {1 − γi ln Ci3 } RW

when ρ < 1 − λi (61)

The friction function, bi (ρ), is discussed in details elsewhere [20], and the empirical model used here is the same as the one used in the single-solute model:
  Ei /RW  (b ) exp when ρ < 1 − λi Faxen i 1−ρ bi (ρ) =  ∞ when ρ ≥ 1 − λi (62) where Ei is the frictional parameter for solute i and (bFaxen )i is given by the Faxen equation [21,25], −1

(bFaxen )i = (1 − 2.104λi + 2.09λ3i − 0.95λ5i )

(63)

Eq. (62) allows the friction function to vary markedly along the radial position inside the pore, and the equation reduces to the Faxen equation when Ei becomes zero, which is the case for very large pores (not practiced in reverse osmosis). Eq. (62) implies that solute diffusivity inside a membrane pore can attain very small values when the parameter Ei is of the same (or larger) order of magnitude as (than) the pore size, RW . 2.1.6. Summary of model assumptions and calculation procedure In this paper, a two-dimensional mathematical model is developed, for transport of multi-solute liquid solutions through RO membranes, based on the following assumptions: • The membrane skin layer is porous, and transport takes place through the fine pores. The membrane transport is similar to the PSCF mechanism. • Strong affinity may exist between the solute species and the membrane polymer. • The fine pores inside the membrane have perfect cylindrical shape, shown in Fig. 1. • Radial equilibrium exists in the pores as discussed above (Eq. (1)). • The membrane system is isothermal. • The liquid solution inside the fine pores is ideal and Newtonian. • van’t Hoff’s equation for osmotic pressure is valid. • The membrane potential function for solute species can be represented by Eqs. (59) and (58), so that the radial and axial effects are separable.

33

• The friction function can be represented by Eq. (62). • The frictional forces inside the fine pores can be modeled by Spiegler’s frictional model [28]. • The effect of axial gradient of the potential function is ignored, in the equation of solute flux through a pore, Eq. (AI.5) in Appendix I, due to an order-of magnitude analysis of all terms of the equation. • The effect of solute–solute frictional forces is ignored, in the equation of solute flux through a pore, Eq. (AI.5) in Appendix I, due to an order-of magnitude analysis of all terms of the equation. The calculation procedure, knowing the membrane parameters (RW , θ i , γ i , Ei ), is as follows: 1. Determine the physicochemical properties of the solvent and solute species needed, such as diffusivities of solute species in the solvent and partial molal volumes of the solute species. 2. Make initial guess for the permeate concentrations, Ci3 , of all solute species. 3. Use Eq. (59) for calculating numerical values of membrane potential function. 4. Use Eq. (62) for calculating numerical values of friction function. 5. Make an initial guess for the adsorbed layer thickness, h. A good initial guess would be an average of all molecular diameters of solute species. For next guess, see step 8. 6. Solve the nonlinear differential equation (35) to find the dimensionless velocity profile in a pore. Other equations, such as Eqs. (19), (24), (25), (27), (28) and (33), are needed to solve Eq. (35). For the numerical technique used in this work, see the end of this section. 7. Calculate the permeate concentrations, Ci3 , of all solute species using Eq. (51). Other equations, such as Eqs. (42), (44) and (45), are needed to solve Eq. (51). 8. Calculate the flux ratio from Eq. (49b). Then use this value to determine a new value for h using Eq. (55). If the difference between the two h values is outside a tolerance, then, modify the h value and return to step 4 and repeat the above steps until convergence is achieved. 9. Compare these calculated Ci3 values of step 7 with those guessed in step 2. Modify your guessed values and repeat steps 3–8 until convergence is achieved between the results. 10. Calculate separation achieved using Eq. (54). By now, both separation and flux ratio are known. Regarding the numerical technique used in this work to solve the nonlinear differential equation (35), previously [29], the method of orthogonal collocation was successfully applied to the governing equations of membrane transport, which proved to give very accurate results. However, recently [30], it has been found that the method of alternatingdirection implicit (ADI) can be used very efficiently with accurate results and CPU time much less than that needed by the orthogonal collocation method. Therefore, the ADI

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

34

Table 1 Physicochemical parameters for each solute Solute j (m3 /kmol)

νj DjB (×109 m2 /s) Rj (×1010 m)

Benzene 0.089 1.096 2.23

Table 2 Optimized parameters based on single-solute experimental data Toluene 0.106 0.968 2.52

Solute j (×1010

θj m) Ej (×1010 m) γj Average RW (×1010 m)

Benzene

Toluene

−29.16 18.93 0.0414

−34.08 35.06 0.1118 12.19

method has been employed in the present research to obtain simulation results. Both numerical techniques have been tested for limiting cases of single-solute under the conditions of zero potential function and zero friction function where analytic solutions exist for the simplified model, and the numerical solutions exactly matched the analytic solutions [20].

which combines the Gauss (Taylor series) method and the method of steepest descent. The unknown model parameters are adjusted to minimize the sum of squares between the model and experimental separation values (f ). Therefore, the objective function to be minimized was defined as the square of differences between model and experimental separation values over all experimental conditions:

3. Results and discussion

Fobjective =

N

  2 (fmodel,k − fexp ’l,k )

(64)

k=1

Simulation study results are presented in this section. Experimental evaluation of the proposed model will be presented in another paper. The solute–membrane affinity systems of single-solute (such as, benzene–water–CA and toluene–water–CA) were studied before at room temperature [9–12]. The benzene and toluene systems have physicochemical properties as shown in Table 1. The experimental data, for toluene–water–CA and benzene–water–CA systems at room temperature, were obtained by Dickson et al. [11,12] at a feed concentration range of 0.4–4.0 mol/m3 for toluene in water, and 0.5–3.5 mol/m3 for benzene in water, and applied pressure range of 690–6900 kPa. The CA membrane samples are cut from flat-sheet membranes and installed in RO cells, and the RO experiments were run after the membrane samples were compacted for several hours to reach steady-state permeation fluxes. The experimental setup (including six RO permeation cells) and procedure are described elsewhere [9–12]. The membrane used was asymmetric cellulose acetate. The experimental data for the above single-solute systems [11,12] were used to optimize the parameters of both singlesolute model [20] and the present multi-solute model for such systems (that is, the present model was tested for the limiting case of single solute). Of course, the present form of potential function was used in the single-solute model. In this way, the validity of the numerical scheme could be tested for the present model. It was found that both models yielded the same optimized values for their parameters. Therefore, the numerical scheme was employed correctly in this study, and the present model correctly reduces to the single-solute model of 1993 when the present potential function, Eq. (59), is employed. The results of the optimized parameters (RW , θ i , Ei and γ i ) are shown in Table 2. A nonlinear parameter estimation code, called UWHAUS [22], was employed to optimize model parameters in this study. The optimizer is based on Marquart’s method [31]

where f is the separation value for the kth experimental conditions, and the counter k runs from first to last (Nth) experimental conditions. The experimental data together with the fitted model for the toluene–water–CA are shown in Fig. 2. The fitted model is based on minimization of the objective function defined by Eq. (64). Also, the experimental data together with the fitted model for the benzene–water–CA are presented in Fig. 3. As it can be seen, in Figs. 2 and 3, the separation will decrease with an increase in applied pressure. It should be noted that the scatter in the experimental data is because of the fact that the feed concentration was varied in the range of 0.4–4.0 mol/m3 (for toluene) and 0.5–3.5 mol/m3 (for benzene) and, secondly, because of experimental error in dealing

Fig. 2. Experimental data [12] of toluene separation vs. operating pressure for aqueous toluene–CA system with toluene concentration in a range of 0.40–4.0 mol/m3 . The solid curve is the best-fit model (with optimized parameters given in Table 2) through the experimental data at an average concentration of 2.0 mol/m3 .

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

Fig. 3. Experimental data [12] of benzene separation vs. operating pressure for aqueous benzene–CA system with toluene concentration in a range of 0.50–3.5 mol/m3 . The solid curve is the best-fit model (with optimized parameters given in Table 2) through the experimental data at an average concentration of 2.0 mol/m3 .

with such low-concentration systems. Also, the experimental data in Figs. 2 and 3 show that as the feed concentration increases, in general, the solute separation will increase at a fixed operating pressure. Next, the above single-solute parameters, in Table 2, were used to simulate mixtures of benzene–toluene in water for the same type of CA membrane under different operating conditions. The performance of the membrane is expressed in term of separation, f , and flux ratio, NT /NP , where NP is pure water flux (i.e., no solute species in the system) through the membrane. The simulation results for the multi-component systems are shown in Figs. 4 and 5, in terms of separation f versus operating pressure, for aqueous mixtures of toluene–benzene

Fig. 4. Simulation results of toluene separation vs. operating pressure for aqueous toluene–benzene–CA system with toluene at 2.0 mol/m3 and different benzene concentrations. Benzene concentration varies from 0.5 to 2.0 mol/m3 . Pure toluene curve is also included.

35

Fig. 5. Simulation results of benzene separation vs. operating pressure for aqueous toluene–benzene–CA system with benzene at 2.0 mol/m3 and different toluene concentrations. Toluene concentration varies from 0.5 to 2.0 mol/m3 . Pure benzene curve is also included.

at low concentrations. The optimized parameters for toluene system (Table 2) are used in the present model to get the results of Fig. 4. The toluene concentration is kept constant at 2.0 mol/m3 , and benzene concentration varies from 0 to 2.0 mol/m3 . As suggested by the model, in general, separation will decrease with increase in applied pressure, and increasing the feed concentration of similar solute species (i.e., benzene) will increase the separation of toluene. Also, the optimized parameters for benzene system (Table 2) are used in the model to get the results in Fig. 5. The benzene concentration is kept constant at 2.0 mol/m3 , and toluene concentration varies from 0 to 2.0 mol/m3 . The same trends of Fig. 4 are seen in Fig. 5. Increasing the feed concentration of similar solute species (i.e., toluene) will increase the separation of benzene. These effects will be discussed more in the following section. Next, the effect of applied pressure and feed concentration on flux ratio are investigated. For the toluene–water–CA system, where experimental data are presented for toluene separation versus pressure in Fig. 2, the experimental data of flux ratio (NT /NP ) versus applied pressure are presented, together with the model prediction values, in Fig. 6. These are the results for a singlesolute (toluene) system, and, in general, the flux ratio will decrease with increase in applied pressure. It should be noted that the main reason for the scatter of experimental data is that these data were obtained for different feed concentrations in a range of 0.40–4.0 mol/m3 . The predicted curve by the model is for a constant feed concentration of 2.0 mol/m3 . It should be noted that the model prediction for flux ratio is based on the objective function in Eq. (64), which considers the separation values only and does not take flux ratio into account for optimization.

36

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

In general, the following simulation results were obtained:

Fig. 6. Experimental data [12] of flux ratio (NT /NP ) vs. operating pressure for aqueous toluene–CA system with toluene concentration in a range of 0.40–4.0 mol/m3 . The dashed curve is the best-fit model (with optimized parameters given in Table 2) through the experimental data at an average concentration of 2.0 mol/m3 .

Next, the simulation results for the effect of feed composition on flux ratio are investigated in Fig. 7. The results of flux ratio (NT /NP ) versus operating pressure for aqueous toluene–benzene–CA system with toluene concentration at 2.0 mol/m3 and different benzene concentrations. Benzene concentration varies from 0.5 to 2.0 mol/m3 . In general, for most of the operating pressure range, increasing the feed concentration will decrease the flux ratio. However, there seems to be another zone in the pressure range where an opposite phenomenon is suggested by the model. This is in a high-pressure zone where the flux ratio may increase by an increase in the feed concentration. This point will be discussed further in the following section.

Fig. 7. Simulation results of flux ratio (NT /NP ) vs. operating pressure for aqueous toluene–benzene–CA system with toluene at 2.0 mol/m3 and different benzene concentrations. Benzene concentration varies from 0.5 to 2.0 mol/m3 . Pure toluene curve is also included.

• Effect of operating pressure towards membrane performance: Increasing operating pressure will decrease the separation of both benzene and toluene, which is a characteristic of the presence of strong solute–membrane affinity for any low feed concentration. This is because an increased convective velocity, at higher pressure, increases the mobility of solute molecules inside pores due to higher convective shear forces against the attraction potential force between the solute and the membrane structure. Generally, the solute–water separation can be positive, negative or zero depending on specific operating conditions of the system. The phenomenon of decreasing separation with increase in applied pressure has been observed by many researchers as stated in the Introduction (see for example [16]). The permeation flux, and flux ratio, will also decrease with increase in operating pressure. This is probably because of a decrease of adsorbed solute layer, h, and consequent convection of the solute molecules along the fluid inside the membrane pores at the higher applied pressures. Again, this phenomenon has been experimentally observed by many researchers (see for example [16]). • Effect of increasing solute concentration towards membrane performance: According to the above simulation results, increasing the concentration of benzene (or toluene), and hence increasing the total concentration, will increase separation of toluene (or benzene). This is because of the weakening of the membrane potential function (as the negative potential function becomes less negative at higher feed concentrations), and hence reduction of the attraction forces between the solute and the membrane structure. Therefore, the presence of benzene will enhance the separation of toluene because increasing feed concentration of benzene with constant toluene concentration (decreasing mole fraction of toluene), which increases the total feed concentration of solutes, will increase the separation due to the weakening of the negative attraction potential field. The separation of solute species (i.e., toluene or benzene) in mixture is higher that in single-solute solution. This effect of increasing separation at higher feed concentrations has been observed by some researchers. For example, Ahmad and Tan [16] studies the effect of pressure on separation of two-solute systems of phenol/pnitrophenol aqueous mixtures and 3-chlorophenol/pnitrophenol aqueous mixtures at different concentrations (0.42–1.48 mol/m3 ) and a range of 500–1500 kPa operating pressure. Their experiments explicitly shows higher separation is achieved at higher feed concentrations. Regarding the effect of increasing the feed concentration of solute species on permeation flux (NT ),

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

or flux ratio (NT /NP ), in general, the flux ratio will decrease with increase in the feed concentration. This typical behavior has been observed by some researchers (e.g. [10–12]). This is most probably due to an increase in the adsorbed layer thickness, h. However, as can be seen in Fig. 7, there could be another zone in the range of applied pressure that the flux ratio may increase with increase in feed concentration. If so, this latter phenomenon will happen at higher pressures. This may be due to increased convective shear forces at the higher pressures against the attraction potential force between the solute molecules and the membrane polymer in the adsorbed layer which, hence, lead to reduction in the thickness of the adsorbed layer at the higher pressures. That is, it is possible that the balance of membrane attraction forces and fluid drag forces becomes upset at higher pressures such that the drag forces prevails the attraction forces, and, therefore, the adsorbed layer becomes thinner and allows for more passage of the fluid flux through the membrane pores. For the time being, the authors could not find any experimental data in the published literature to support this possible phenomenon, and this could be investigated in a well-designed experimental research in future.

Appendix A. Force balance on the solute in a membrane pore A force balance on solute i inside the pore yields       n   Fi (r, z) = − FiB (r, z)+FiM (r, z) + Fij (r, z)       j=1 j = 1

(AI.1) where FiB , FiM , and Fij are frictional force between the solute i and solvent, membrane, and solute j, respectively, which are given by the frictional model [28] as FiB (r, z) = −χiB [ui (r, z) − uB (r)] FiM (r, z) = −χiM [ui (r, z) − 0] = −χiM

Acknowledgement The authors wish to acknowledge the research grant provided by Intel Technology Sdn. Bhd. (Malaysia) for providing financial support for this research.

(AI.2) Ji,z (r) Ci (r, z)

Fij (r, z) = −χij [ui (r, z) − uj (r, z)]

(AI.3) (AI.4)

where χij is the frictional constant between components i and j. Combining Eqs. (AI.1)–(AI.4) with Eq. (15) for the driving force Fi (r, z) gives  πi (ρ, ξ) ∂P(ρ, ξ) 1 DiB  ∂πi (ρ, ξ) Ji,z (ρ) = − − − υi τRT bi (ρ)  ∂ξ RT ∂ξ

4. Conclusions A mathematical model has been developed for multicomponent reverse osmosis systems in which strong affinity exists between solutes and membrane polymer. The model parameters optimized based on single-solute systems (such as toluene or benzene) were used to simulate multi-component system of aqueous benzene–toluene–CA system at different operating conditions. The simulation studies suggest that both benzene and toluene separation will decrease with increase in operating pressure, with toluene possessing higher separations, and will increase with increase in feed concentration. Also, permeation flux (or flux ratio) will decrease with increase in operating pressure and will decrease, in general, with increase in solutes concentration of feed depending on the operating pressure. An opposite phenomenon of increasing flux ratio with increase in feed concentration could happen, according to the present model, at high-enough pressures presumably due to a decrease in thickness of adsorbed layer of solutes at the higher pressures.

37

¯ D ∂Φi (ρ, ξ) − πi (ρ, ξ) DiB ∂ξ  n  χij + πi (ρ, ξ) τJj,z (ρ)  πj (ρ, ξ) + α(ρ)πi (ρ, ξ)

(AI.5)

(j=1)=1

where ρ and ξ are the dimensionless coordinates defined by Eqs. (7) and (8), Φi (ρ, ξ) is defined by Eq. (6), α(ρ) is the dimensionless fluid velocity defined by Eq. (18), bi (ρ) is the friction parameter defined by Eq. (17). In deriving Eq. (AI.5), the following relationship [28]: DiB =

RT χiB

(AI.6)

and the equation for solute velocity: ui (r, z) =

Ji,z (r) Ci (r, z)

(AI.7)

have been used. Eq. (AI.5) can be further simplified, by ignoring the effect of the potential gradient and solute–solute frictional forces, to the following form:  1 DiB ∂πi (ρ, ξ) πi (ρ, ξ) ∂P(ρ, ξ) Ji,z (ρ) = − − − υi τRT bi (ρ) ∂ξ RT ∂ξ ¯ D + α(ρ)πi (ρ, ξ) DiB

 (AI.8)

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

38

which is given as Eq. (16). The justification for the above simplifying assumptions will be discussed later. The effect of axial gradient of the potential function is ignored in the equation of solute flux through a pore, Eq. (AI.5), due to an order-of magnitude analysis of all terms of the equation. This can be justified by employing the numerical values of model parameters in determining the value of the potential function and its gradient in axial direction. Preliminary calculations reveal that this term is about 10 times smaller than the remaining terms. Also, the effect of solute–solute frictional forces is ignored in the equation of solute flux through a pore, Eq. (AI.5), due to an order-of magnitude analysis of all terms of the equation. This is true especially that the term includes τ (the effective length of a pore) and the solute flux itself which have small numerical values. The χij value can be determined by using an equation similar to Eq. (AI.6): χij = RT/Dij . Preliminary calculations reveal that this term is about 100 times smaller than the remaining terms.

Appendix B. Derivation of velocity profile The force balance on the elemental fluid of Fig. 1 consists of the following parts: (i) The net force due to difference in pressure  
  ∂P  (2πr dr) (P|z − ( P|z + dz ∂z z = 2πr dr dz

∂P(r, z) ∂z

(AII.1)

(ii) The net force due to viscous shear stresses, using Newton’s law of viscosity,  ∂2 uB (r) (2πr dz){(P|r − (P|r+dr } = η 2πr dr dz ∂r 2  duB (r) + (2πr dr dz) (AII.2) dr (iii) The net force due to the friction force between solutes and the pore wall, n

FiM (2πr dr dz)Cj (r, z) = (2πr dr dz)

j=1

×

n

[−χjM (r)uj (r, z)]Cj (r, z) = −(2πr dr dz)

j=1

×

n

χjM (r)Jj,z (r)

(AII.3)

becomes: n

FjM (ρ, ξ)(2πρ dρ dξ)Ci (ρ, ξ) = −(2πρ dρ dξ)

j=1

×

j=1

 χjM (ρ)

j=1

 ×

DiB [mi α(ρ) + ωi (ρ)] bi (ρ) τRT

πi2 − K∗ (ρ)πi3 πi2 + [m α(ρ)+ωi (ρ)] i e i −1





e−Φi (ρ,0)

(AII.4)

Then the differential equation for the fluid velocity can be obtained by adding up all the three contributions of the force balance as given by Eq. (35).

Nomenclature bj (ρ) friction function defined by Eq. (17) C molar density of solution (kmol/m3 ) Cj (ρ, ξ) concentration of solute j inside a pore at position ρ and ξ (kmol/m3 ) dj parameter for solute j defined by Eq. (52) ¯ D average diffusivity (of solutes in free solution) defined by Eq. (19) (m2 /s) DjB diffusivity of solute j in free solvent (m2 /s) Ej friction parameter for solute j in Eq. (62) (m) fj theoretical separation of solute j, defined by Eq. (54) Fobjective objective function defined by Eq. (64) J flux through a single pore (kmol/m2 s) ∗ Kj (ρ) ratio of local partition coefficients at the ends of a pore for solute j, defined by Eq. (33) l1 definite integral defined by Eq. (42) l2i definite integral defined by Eq. (44) (kPa) mi parameter defined by Eq. (27) n number of solute species in feed solution N number of experimental observations in Eq. (64) NB solvent (water) flux through membrane (kmol/m2 s) NP pure solvent (water) flux through membrane (kmol/m2 s) NS total flux of all solutes through membrane (kmol/m2 s) NT total permeation flux through membrane (kmol/m2 s) P pressure difference across membrane (kPa) R gas constant (kJ/kmol K) Rj average radius of solute j (m) T temperature (K)

j=1

where χiM (r) is defined in Eq. (AI.3). Substituting for Ji,z (r) from Eq. (34) and substituting for dimensionless axial and radial coordinates, Eq. (AII.3)

Greek letters α(ρ) fluid velocity inside a pore, defined by Eq. (18) β1 parameter defined by Eq. (36)

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

second potential parameter for solute i in Eqs. (57) and (59) ∆ vector of gradient in Eq. (20) ε fractional pore area of membrane η solution viscosity (kPa s) θi first potential parameter for solute i in Eqs. (57) and (59) (m) θ1 potential parameter in Eq. (56) (m) λj ratio of solute j molecular radius to pore radius νj partial molal volume of solute j (m3 /kmol) ξ axial coordinate inside a pore πj2 contribution to osmotic pressure of solution by solute j at feed–membrane interface (kPa) πj3 contribution to osmotic pressure of solution by solute j at permeate–membrane interface (kPa) π2 osmotic pressure of solution at feed– membrane interface (kPa) ρ radial coordinate inside a pore σ j2 (ρ) local Staverman (reflection) coefficient by solute j at feed–membrane interface defined by Eq. (24) σ j3 (ρ) local Staverman (or reflection) coefficient by solute j at permeate–membrane interface defined by Eq. (25) τ average pore length taking tortuosity into account (m) Φi (ρ, ξ) potential function for solute i, dimensionless ωi (ρ) parameter defined by Eq. (28) γi

Subscripts B solvent exp’l experimental i (or j) solute species i (or j) iM solute i in membrane pore jB solute j in solvent k counter of experimental conditions in Eq. (64) model by model P pure solvent (pure water) S solutes T total 2 feed at membrane interface 3 permeate

References [1] H.K. Lonsdale, U. Merten, M. Tagami, Phenol transport in cellulose acetate membranes, J. Appl. Polym. Sci. 11 (1967) 1807. [2] T. Matsuura, S. Sourirajan, Physicochemical criteria for reverse osmosis separation of alcohols, phenols and monocarboxylic acids in aqueous solutions using cellulose acetate membranes, J. Appl. Polym. Sci. 15 (1971) 2905. [3] J.E. Anderson, S.J. Hoffman, C.R. Peters, Factors influencing reverse osmosis rejection of organic solutes from aqueous solutions, J. Phys. Chem. 76 (1972) 4006.

39

[4] T. Matsuura, S. Sourirajan, Reverse osmosis separation of phenols in aqueous solutions using cellulose acetate membranes, J. Appl. Polym. Sci. 16 (1972) 2531. [5] T. Matsuura, S. Sourirajan, M.E. Beduas, Reverse osmosis separation of single and mixed alcohols in aqueous solutions using porous cellulose acetate membranes, J. Appl. Polym. Sci. 18 (1974) 567. [6] T. Matsuura, P. Blais, J.M. Dickson, S. Sourirajan, Reverse osmosis separations for some alcohols and phenols in aqueous solutions using aromatic polyamide membranes, J. Appl. Polym. Sci. 18 (1974) 3671. [7] J.M. Dickson, T. Matsuura, S. Sourirajan, Transport characteristics in the reverse osmosis system p-chloropenol–water–cellulose acetate membrane, Ind. Eng. Chem. Process. Des. Dev. 18 (4) (1979) 641. [8] H.-G. Burghoff, K.L. Lee, W. Pusch, Characterization of transport across cellulose acetate membranes in the presence of strong solute–membrane interactions, J. Appl. Polym. Sci. 25 (1980) 323. [9] D.R. Lloyd, M. Babai-Pirouz, J.M. Dickson, Cyclic hydrocarbon– water separations using a pressure-driven membrane separation process, J. Sep. Technol. 3 (4) (1982) 21. [10] J.M. Dickson, M. Babai-Pirouz, D.R. Lloyd, Aromatic hydrocarbon– water separations by a pressure-driven membrane separation process, Ind. Eng. Chem. Process. Des. Dev. 22 (1983) 625. [11] J.M. Dickson, D.R. Lloyd, The effect of solute–membrane affinity on cyclic hydrocarbon–water transport in pressure-driven membrane separation processes, J. Membr. Sci. 21 (1984) 33. [12] J.M. Dickson, Reverse osmosis in the presence of strong solute–membrane affinity, Ph.D. Thesis, Virginia PolyTechnic Institute and State University, 1985. [13] P.J. Connell, J.M. Dickson, Modeling of reverse osmosis separations with strong solute–membrane affinity at different temperatures using the finely porous model, J. Appl. Polym. Sci. 35 (5) (1988) 1129. [14] Y. Fang, S. Sourirajan, T. Matsuura, Reverse osmosis separation of binary organic mixtures using cellulose acetate butyrate and aromatic polyamide membranes, J. Appl. Polym. Sci. 44 (1992) 1959. [15] V.V. Goncharuk, D.D. Kucheruk, V.M. Kochkodan, V.P. Badekha, Removal of organic substances from aqueous solutions by reagent enhanced reverse osmosis, Desalination 143 (2002) 45. [16] A.L. Ahmad, K.Y. Tan, Reverse osmosis of binary organic solute mixtures in the presence of strong solute–membrane affinity, Desalination 165 (2004) 193. [17] M. Soltanieh, W.N. Gill, Review of reverse osmosis membranes and transport models, Chem. Eng. Commun. 12 (1981) 279. [18] H.K. Lonsdale, U. Merten, R.L. Riley, Transport properties of cellulose acetate osmotic membranes, J. Appl. Polym. Sci. 9 (1965) 1341. [19] S. Sourirajan, Reverse Osmosis, Academic Press, New York, 1970. [20] H. Mehdizadeh, J.M. Dickson, Modeling of reverse osmosis in the presence of strong solute–membrane affinity, AIChE J. 39 (3) (1993) 434. [21] H. Mehdizadeh, J.M. Dickson, Theoretical modification of the surface force-pore flow model for reverse osmosis transport, J. Membr. Sci. 42 (1989) 119. [22] H. Mehdizadeh, J.M. Dickson, Evaluation of the surface force-pore flow and modified surface force-pore flow models of reverse osmosis transport, Chem. Eng. Commun. 103 (1991) 65. [23] S. Jain, K. Gupta, Analysis of modified surface force-pore flow model with concentration polarization and comparison with Spiegler–Kedem model in reverse osmosis systems, J. Membr. Sci. 232 (2004) 45. [24] G. Jonsson, Methods for determining the selectivity for reverse osmosis membranes, Desalination 24 (1978) 19. [25] J.L. Anderson, D.M. Malone, Mechanism of osmotic flow in porous membranes, Biophys. J. 14 (1974) 957. [26] A. Higuichi, T. Nakagawa, Membrane potential and ion transport in inhomogeneous ion-exchange membranes, J. Chem. Soc., Faraday Trans. 85 (1989) 3609.

40

H. Mehdizadeh et al. / Journal of Membrane Science 267 (2005) 27–40

[27] S.W. Thiel, D.R. Lloyd, J.M. Dickson, Physicochemical interpretation of the behavior of a pressure-driven membrane separation process, ACS Symp. Ser. 120 (1985). [28] K.S. Spiegler, Transport processes in ionic membranes, Trans. Faraday Soc. 54 (1958) 1408. [29] H. Mehdizadeh, J.M. Dickson, Solving nonlinear differential equations of membrane transport by orthogonal collocation, Comput. Chem. Eng. 14 (2) (1990) 157.

[30] H. Mehdizadeh, Y.C. Chong, Numerical simulation of velocity and concentration profiles in reverse osmosis membranes using the numerical methods of orthogonal collocation, successive overrelaxation, and alternating-direction implicit, Eng. J. Univ. Qatar 17 (2004) 51. [31] D.L. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, J. Soc. Ind. Appl. Math. 2 (1963) 431.