A mathematical study of the evolution of fouling and operating parameters throughout membrane sheets comprising spiral wound modules

Chemical Engineering Journal 187 (2012) 222–231

Contents lists available at SciVerse ScienceDirect

Chemical Engineering Journal journal homepage: www.elsevier.com/locate/cej

A mathematical study of the evolution of fouling and operating parameters throughout membrane sheets comprising spiral wound modules Margaritis Kostoglou a , Anastasios J. Karabelas b,∗ a b

Division of Chemical Technology, Department of Chemistry, Aristotle University, Greece Chemical Process Engineering Research Institute, Centre for Research and Technology – Hellas, 6th km Charilaou-Thermi Road, GR 570 01, Thermi-Thessaloniki, Greece

a r t i c l e

i n f o

Article history: Received 24 October 2011 Received in revised form 29 December 2011 Accepted 3 January 2012 Keywords: Fouling spatial and temporal evolution on membrane sheets Spiral wound membrane modules Theoretical analysis Simulation algorithms

a b s t r a c t Local experimental data on key process parameters, throughout membrane sheets of spiral-wound membrane (SWM) elements, are unavailable and very difficult to obtain; therefore, theoretical models are essential to better understand the evolution of fouling, which is helpful for improving design and operation of these modules. The mathematical model, formulated in a fundamental manner herein, takes into account the spatial flow field non-uniformities inside narrow channels with spacers (simulating SWM elements) and their interaction with the evolving deposits. A realistic meso-scale transient model, typical of the early stages of organic membrane fouling, is employed. Appropriate non-dimensionalization of model equations allows to treat separately the two (practically significant) cases of constant inlet pressure and constant average permeate flux. Several simplified cases are examined, for specific dimensionless parameter values, and asymptotic solutions are derived. The results of this analysis have guided the development of very efficient numerical techniques, most useful for future extensions towards the development of a comprehensive modeling tool, appropriate for module design calculations and for assessing detailed experimental data. © 2012 Elsevier B.V. All rights reserved.

1. Introduction In modern desalination and water treatment plants, the spiral wound membrane (SWM) element is the most significant component, where the water filtration process takes place [1]. Six to eight SWM elements are usually connected in series within a pressure vessel, so that process conditions significantly vary along this vessel as well as throughout the individual membrane sheets (comprising a SWM) of active surface area of order 1 m2 . Key process parameters, varying as a result of flow through the narrow SWM channels and the membrane, include local permeate flux, trans-membrane pressure (TMP), retentate flow velocity. Fouling is considered as the main problem that degrades SWM element performance, negatively affecting energy losses, permeate quality and overall process efficiency. The development of a fouling layer on membrane sheets is determined by the aforementioned locally varying flow parameters; it is, therefore, characterized by spatial non-uniformity. To improve SWM element design and operation, adequate knowledge (not available at present) is necessary on the temporal and spatial variation of the above parameters across pressure vessels.

∗ Corresponding author. Fax: +30 2310 498189. E-mail address: [email protected] (A.J. Karabelas). 1385-8947/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cej.2012.01.028

Problems encountered in obtaining an improved understanding of the transport processes and fouling development inside SWM elements can be briefly summarized as follows. At the experimental level, two main problems are identified: (a) The test sections in laboratory studies are quite small (usually employing a membrane area two orders of magnitude smaller than a typical membrane sheet), thus simulating only “local” conditions of the aforementioned main parameters. (b) The compact design of SWM elements (with channels of very narrow gap) has not permitted so far any kind of local measurements throughout a membrane sheet that would shed light on prevailing conditions inside such modules. At the theoretical level, progress is hindered by the complexity of flow field geometry and of the related transport phenomena; it is noted in particular that: (a) Models of the flow field at steady state (with respect to the deposit evolution time scale) are available providing some useful basic information; however, modeling fouling rates is much less developed, as significant uncertainties exist on prevailing mechanisms. (b) Models including the temporal evolution of fouling, throughout a membrane sheet, are not available. A major difficulty holding

M. Kostoglou, A.J. Karabelas / Chemical Engineering Journal 187 (2012) 222–231

back progress in this area lies in the fact that, even for relatively simple types of fouling, quantitative relationships linking fouling rates with the main flow parameters (local TMP, flux, cross-flow velocity) are not available. For these reasons, the problem of predicting evolution of fouling rates and of concomitant flux decline, throughout SWM sheets, is essentially untouched to this day. In fact, a comprehensive review on SWM elements [2] concludes that “. . .modeling of flux decline over time due to depositions on the membrane surface remains one of the most intractable issues for membrane processes”. This study aspires to contribute towards improving this situation, taking advantage of recent relevant theoretical and experimental work by the authors and their collaborators [3–8]. The scope of the present paper is to mathematically model in a fundamental manner the evolution of deposit on a sheet, typical of SWM module, by taking into account the spatial non-uniformities of the flow field and their interaction with the evolving deposit. Several types of membrane fouling occur in practical water treatment operations; i.e., scaling due to sparingly soluble salts [9,10], organic gel-layer formation [7,11–13], biofouling [14], colloidal particle deposition [15,16]. In general, the effect of the deposit is to reduce the effective permeability of both the membrane and the narrow retentate channel, with obvious negative consequences. The relative significance of this effect on channel and membrane permeability depends on the type and extent of fouling. The present work is focused on relatively thin fouling deposits growing on an initially clean membrane. This is a common (and relatively simple from the modeling point of view) mode of fouling, representative of organic layer formation [7,11]; moreover, it is known that this type of organic fouling is due to growth of a compressible gel layer (also called “cake”) on the membrane [6,7]. There are several theoretical studies treating this particular mode of fouling essentially as a dead-end filtration process, which is a good assumption for the case of insignificant pressure drop in the retentate and permeate sides of the SWM. However, during the practical operation of a SWM module, the pressure drop at the retentate and permeate sides promotes a 2-dimensional distribution of the trans-membrane pressure, thus creating a distribution of the wall permeate flux throughout the membrane sheet [2,17]. In turn, this leads to a distribution of the fouling layer thickness and correspondingly of the local membrane permeability. Moreover, the deposit non-uniformity interacts with the mesoscopic flow field in the SWM which must be resolved in order to predict the temporal evolution of the deposit. Existing, fundamentally formulated, models for membrane module operation are capable of representing in detail specific aspects of the process (e.g. the flow details in narrow channels with spacers [3]) and can be characterized as ‘learning’ models [18]. These models cannot reproduce the full module operation so they cannot be used as design, optimization and control tools. However, typical mesoscopic models [2,5,17] are considered appropriate for this purpose, based on mass and momentum balances as well as on constitutive laws extracted from experimental results or from detailed subscale calculations. Simplified analytical solutions of this type of models are available in the membrane literature (e.g. [19]). It is pointed out that such solutions usually facilitate understanding of the physical processes studied, as suggested in [20], and in many cases guide the development of problem specific numerical techniques. In view of its advantages, the above approach is taken in this study. It should be pointed out here that the scope of the present work is not to develop models of the actual desalination process but efficient tools for interpretation of experimental data from well-controlled, idealized experiments focused on specific phenomena (organic layer growth in this particular case) encountered in actual desalination processes. In this context, numerical techniques for the treatment

223

of the phenomena under study are developed for inclusion as a module in a generalized desalination simulation tool. The structure of the present paper is as follows: first, a mesoscopic transient model for deposit evolution on a membrane sheet is developed and properly non-dimensionalized. Then several simplified cases for specific values of parameters are examined and asymptotic solutions are derived. Exploiting the results of this analysis, very efficient numerical techniques are developed. Several indicative results are presented in the last section. 2. Problem formulation 2.1. Retentate channel For the purpose of this study, a flat geometry of the intermembrane flow channels is considered, as shown in Fig. 1. Although the picture outlined in the preceding section seems rather simple, the complicated geometric structure of the SWM (in particular at the retentate/concentrate side [21,22]) renders the modeling of the flow field difficult. A direct handling of the problem with simultaneous consideration of all the spatial and temporal scales is not possible, so scale separation procedures (employed here) are helpful. The small scale closure relation is the “constitutive” expression relating local pressure drop in the retentate channel to the local cross-sectionally averaged flow velocity. Such a constitutive relation can be obtained either by fitting appropriate expressions to experimental pressure drop vs. flow rate data, or from subscale CFD simulations. In the latter case several approaches have been taken [23–25], with prevalent the 3-dimensional unit cell one [3]. Since this work is focused on the deposit evolution, the fouling model is outlined first. A rather simple, yet realistic, cake filtration model is considered, whereby the dissolved organic material is transferred with the wall velocity to the membrane surface, sticking onto the developing fouling layer; all the organic matter corresponding to the permeate is considered to deposit on the membrane. Rather extensive experimental work in this Laboratory [6,7] and elsewhere (e.g. [12,13]) shows that this type of deposit layer (cake) formation mechanism adequately describes RO/NF membrane fouling by common organic matter, such as polysaccharides and their mixtures with humic acids; indeed, in the initial stage of development which is of practical interest [7], such deposits are thin coherent, gel-type, layers formed by the organic matter present in the fluid

Fig. 1. Schematic of the membrane element geometry (a) flow cross section (b) top view.

224

M. Kostoglou, A.J. Karabelas / Chemical Engineering Journal 187 (2012) 222–231

permeating the membrane. No re-entrainment due to cross flow and cake enhanced concentration polarization [26] are considered. In addition, the realistic assumption is made that the gel layer is very thin, compared to the retentate channel gap, so that the effective permeability of this channel is not altered by the deposit. According to the above, the evolution of the local mass surface density m of the deposit is simply described by the equation dm = uw c dt

(1)

where uw is the normal wall velocity (permeate flux) and c is the concentration of the depositing material present in the liquid. According to the assumptions made above, c is uniform throughout the channel and equal to its value in the feed stream. The relation between uw and trans-membrane pressure p is commonly given as [27]: uw =

1 p (Rm + Rc )

(2)

where  is the liquid viscosity, Rm the clean membrane resistance and Rc is the fouling cake resistance related to the deposit mass through Rc = ˛m, where ˛ is the specific cake resistance. For compressible cakes (as is the case of gel layers examined here), experiments show that ˛ is proportional to (ıp)n where ıp is the difference between the applied pressure (retentate pressure) and a reference pressure and n is the compressibility exponent, i.e. ˛ = ˛O (ıp)n . The increase of the cake resistance due to compressibility effects is due to the fact that a part of the trans-membrane pressure is transformed to cake stresses, so the structure of the cake changes towards smaller porosity. Several models exist in the literature to account for the structural variations across the cake thickness (e.g. [28]). The error of using the pressure p without consideration of solid stresses to relate with ˛ is discussed in [29]. It is important to point out here that, as shown recently [6,7], experimental data on flow resistance evolution for organic fouling layers can be described by considering that ˛ is a function of the pressure drop across the fouling layer, i.e. ˛ = ˛O (p − Rm uw )n [6,7]. This expression implies that the specific resistance tends to increase (starting from a zero value) as the thickness of the cake increases. Although an even more general behavior may be encountered in practice (e.g. the exponent n may be function of cake thickness), the above expression is considered a reasonable approach to a kind of constitutive law for the cake resistance. In fact, this is particularly true for the thin gel layers [6,7] of practical interest in membrane desalination processes. Extensive unit cell CFD simulations have been carried out [3], showing that the pressure drop is related to the superficial, crosssectionally averaged, flow velocity U as follows U2 p = −f1 Re−f2 D L

(3)

where the Reynolds number is defined as Re = UD/ and ,  are the fluid density and viscosity, respectively. The retentate channel pressure is p and the spacer filament diameter is D. The common (net-type) spacers dealt with are comprised of two layers of parallel unwoven filaments; thus, the channel gap is twice the filament diameter. Detailed geometric description of spacers is provided elsewhere [3]. The coefficients f1 and f2 depend on the spacer geometry (pitch to diameter ratio H/D and filament crossing angle ˇ). Generalizing the relation in Eq. (3) by considering a differential pressure drop in an arbitrary direction and substituting for Re leads to  =− U

D1+f2 f1

1−f2 f2 U 1−f2

∇p

(4)

By taking the magnitude of both sides, of the above vectorial identity, solving for U and replacing again in Eq. (4), one obtains



 =− U

1/(2−f2 )

D1+f2 f1 1−f2 f2

 (f2 −1)/(2−f2 ) ∇ p ∇p

(5)

where the velocity vector is averaged along a direction normal to membrane surface. 2.2. Permeate channel For the permeate channel, where a porous material fills the gap, a linear resistance law based on the smaller flow rates and smaller pore characteristic sizes can be typically assumed. The absence of fouling in this side implies a uniform permeability leading to the following relation for the cross sectional average velocity u: =− u

k ∇P 

(6)

where P is the permeate side pressure and k the corresponding permeability. The permeate velocity vector is averaged along the direction normal to the membrane. The differential mass balances, in x and y directions, for the retentate and permeate channels lead to the following equations:  = −uw Retentate side Lrz ∇ U

(7)

 = uw Permeate side Lpz ∇ u

(8)

Combining Eqs. (1)–(8) a system of two coupled partial differential equations for the pressure fields p and P is obtained:



(f2 −1)/(2−f2 )

∇ ∇ p ∇2P = uw =



∇p =

−1/(2−f2 )

D1+f2 f1 1−f2 f2

uw kLpz

uw Lrz

(9) (10a)

p−P (Rm + ˛o (p − P − uw Rm )n m)

(10b)

It is noted that the last equation is implicit with respect to uw . This system needs appropriate boundary conditions. In Fig. 1, the no penetration (zero velocity boundary) conditions for the two boundaries of the retentate side (at y = 0 and y = Ly ) and the three boundaries of the permeate side (at x = 0, x = Lx and y = Ly ) are evident. The pressure at the outflow edge of the permeate side (y = 0) can be set to zero; it is noted that this is not restrictive at all since only pressure differences appear in the equations. The physically appropriate choice for the remaining boundary conditions is to define the retentate channel pressure at x = 0 (pin ) and x = Lx (po ); i.e. to consider the pressure drop as an input parameter and the flow rate as the outcome of the problem solution. In general, the inlet pressure is a function of time and its initial value is denoted as pin,o . Considering that the inlet flow is towards the module and the exiting retentate flows in the same direction, leads to the following expressions for the inlet and outlet velocity profiles:



Uin (y) =

−

 Uout (y) =

D1+f2



f1 1−f2 f2

−

D1+f2 f1 1−f2 f2

∂p ∂x





∂p ∂x

1/(2−f2 )

(11) x=0



1/(2−f2 ) (12)

x=Lx

The average inlet and outlet velocities Uin,av , Uout,av can be Ly obtained by applying the averaging operator 1/Ly 0 . . . dy to the above functions. The liquid recovery fraction R (i.e. ratio of permeate outlet over retentate inlet flow) can be computed as

M. Kostoglou, A.J. Karabelas / Chemical Engineering Journal 187 (2012) 222–231

R = 1 − [Uout,av /Uin,av ]. An alternative way is to use the flux through the membrane, i.e.

Lx Ly R= 0

t¯ =

˛r uwo ct Rm

(22)

with boundary conditions uw dxdy (Ly Lrz Uin,av )

(13)

0

Unlike the clean membrane steady state problem [8], the dynamic fouling problem requires a strategy to deal with variations of operational parameters as the membrane resistance evolves in time. In the present work, the two (practically realistic) cases will be considered; i.e., of fixing po , Uin,av and alternatively Uin,av , R, corresponding to well known constant pressure and constant flux cases of dead-end filtration, respectively. In order to facilitate analytical treatment of the problem, a simplified form of the specific cake resistance ˛ versus pressure dependence is considered; i.e. ˛ = ˛r (1 + (p − P − pin,o )/pin,o ). This expression may be viewed as a generalization of the Taylor expansion of the experimentally obtained form around the exponent n = 0 and it is better applicable to slightly compressible cakes, which is the case for relatively thin fouling layers treated here [6,7]. In this form the effect of compressibility is introduced through the parameter . Moreover, this (somewhat simplified) expression is considered to better represent the behavior of a more general specific resistance ˛ versus effective pressure dependence, while facilitating the mathematical treatment.

∂p¯ = 0 at ∂y¯

y¯ = 0,

∂P¯ = 0 at ∂y¯

y¯ = L¯ y

∂P¯ = 0 at ∂x¯

x¯ = 0,

P¯ = 0 at

The above mathematical problem contains many physical parameters, which must be grouped into a smaller number of dimensionless parameters, convenient for analyzing the equations. All pressures are normalized by pin,o and all lengths and spatial coordinates by Lx . The inlet and outlet velocities are normalized by the velocity Uo defined as Uo =

pin,o D1+f2 f1 1−f2 f2 Lx

1/(2−f2 )

(15)

1

(16)

¯ 1 + (1 + (p¯ − P¯ − 1))m 1

(17)

¯ 1 + (1 + (p¯ − P¯ − 1))m

(18)

where

(27)

x¯ = 1

(28)

The evolution of p¯ in and p¯ o is determined by several constraints which are subsequently discussed. Also a reference value for the flow recovery fraction R can be defined as Ro = uwo Lx /(Lrz Uo ). The above normalization leads to the following expressions for the output variables:



¯ = U¯ in (y)

−

∂p¯ ∂x¯

 ¯ = U¯ out (y)

−

U¯ out,av

1 = L¯ y

B=

1+ε

(29) x¯ =0

∂p¯ ∂x¯

1+ε (30) x¯ =1

L¯ y U¯ out dy¯

(31a)

U¯ in dy¯

(31b)

0



L¯ y

U¯ in,ave =

1 L¯ y

u¯ w,av

1 = ¯Ly

1 ¯ 1 + (1 + (p¯ − P¯ − 1))m

(31c)

1 L¯ y uw dydx

(31d)

0

Ro u¯ w,av U¯ out,av = R=1− ¯ Uin,av U¯ in,av

(31e)

Despite its complicated structure, the problem is determined by only four dimensionless parameters, i.e. A, B, ε,  which have specific physical meaning. The exponent ε denotes the deviation from the steady laminar flow (ε = 0),  reflects the cake compressibility, B is related to the pressure drop in the permeate side and A is related to the recovery fraction. Recent reviews on solution of the governing equations for the clean membrane and some closed form solutions can be found in [5] (two dimensions, only hydrodynamic equations) and [30] (one dimension including solute mass balances). 3. Analysis of typical cases

f2 − 1 ε= 2 − f2 A=

(26)

x¯ = 0

0

¯ p¯ − P¯ dm = dt¯ ¯ 1 + (1 + (p¯ − P¯ − 1))m



(25)

p¯ = p¯ o (t¯ ) at

¯ u¯ w = (p¯ − P)

ε

2 ¯ ∇¯ P¯ = B(P¯ − p)

x¯ = 1

0

¯ = The dimensionless deposit mass density is given as m ˛r m/Rm . It is noted that the specific cake permeability is incorporated into the dimensionless mass which expresses directly the effect of fouling on membrane overall permeability. The dimensionless form of equations for the pressure fields and deposit mass evolution is



(24)

(14)

1 p = Rm  in,o

¯ ∇¯ ∇ p¯  ∇ p¯ = A(p¯ − P)

(23)

p¯ = p¯ in (t¯ ) at

and the permeate flux through the membrane is normalized by the reference value uwo defined as uwo

y¯ = L¯ y

y¯ = 0

2.3. Non-dimensionalization



225

f1 1−f2 f2 D1+f2

Lx2 kRm Lpz

(19)

1/(2−f2 ) 

Lx pin,o

(f2 −1)/(2−f2 )

uwo Lx2 pin,o Lrz

(20)

(21)

The main interest of the present work is in the spatial evolution of the uw and m and their effect on the operating variables of the SWM element. Two modes of membrane operation are considered. The first case is the constant pressure/constant inlet flow mode where p¯ in (t¯ ) = 1, U¯ in,ave = constant. In this case, po results as an output of the solution procedure. The second case is the constant permeate flux/constant inlet flow mode, where u¯ w,av = constant,

226

M. Kostoglou, A.J. Karabelas / Chemical Engineering Journal 187 (2012) 222–231

U¯ in,ave = constant. In this case not only the outlet pressure but the inlet pressure evolution p¯ in (t¯ ) is the outcome of the solution procedure. The direct numerical solution of the above problem is cumbersome. The steady state coupled system of non-linear partial differential equations is stiff due to the non-linearity of the channel flow resistance law. This pseudo steady state problem must be solved repeatedly in time to account for the evolution of the deposit shape (membrane permeability reduction pattern). Imposing a constant inlet flow rate requires an iterative procedure to find p¯ o (t¯ ). A further complication arises for the case of constant wall flux strategy where an additional iteration cycle at each time step is required in order to find the inlet pressure that keeps the average wall flux constant. Fortunately, the mathematical problem can be considerably simplified not only for idealized laboratory conditions but also for conditions met in practical applications of SWM. In the following, several simplified versions of the above set of equations and the corresponding solution techniques will be presented. Subsequently, only dimensionless variables appear so the over-bar is omitted for clarity of presentation.

the deposit grows faster, the permeate flux decreases leading to a reduction of deposition rate. So, there is a tendency of the growth process to eliminate cake thickness non-uniformity imposed by the initial permeate flux non-uniformity. Surprisingly, according to the above analytical solutions, the initial non-uniformity is not completely eliminated but is reduced asymptotically to roughly half of its initial value; i.e. the fractional difference of the deposit thickness between x = 0 and x = 1 varies from  to /2. Constant flux

⎛

m = −1 + ⎝1 + 2 ⎝

m = −1 + (1 + 2t)1/2

(32a)

uw = (1 + 2t)−1/2

(32b)

Constant flux m=t pin

(33a)

1 + t − t = 1 − t

(33b)

(ii) For finite pressure drop along the channel, it can be shown that the pressure profile is linear for any flow resistance exponent. The pressure profile for constant inlet velocity can be described as p = pin − x where  is the constant pressure slope (pin − po ). The governing equations are: dm = uw , dt

uw =

pin − x 1 + (1 + (pin − x − 1))m

(34)

⎞⎞1/2 pin (t  )dt  − xt ⎠⎠

(36a)

0

uw =

t

(1 + 2(

0

pin (t  )dt  − xt))

(36b)

1/2

The function pin (t) must be determined in this case from the solution of the following integral equation

1

(pin (t) − x)

t

Case 1a. Insignificant fluid recovery (A  1) (i) In the particular case of small (fractional) difference between the inlet and outlet pressure there is no driving force for development of any kind of non-uniformity and the problem degenerates to that of dead end filtration. Solutions are given for reference: Constant pressure

t

(pin (t) − x)

3.1. CASE 1. No pressure drop in the permeate side (B  1) In this case P = 0 everywhere, Eq. (17) can be neglected and the problem is reduced to one dimension, since variations appear only in the x-direction.

⎛

0

(1 + 2(

0

pin

(t  )dt 

− xt))

1/2

dx =

1− 2

(37)

Although the solution of the above equation is not easier than that of the original problem (Eq. (34)), it is informative in that Eq. (37) leads to the following value of the asymptotic deposit t non-uniformity between x = 0 and x = 1: t/(2 0 pin (t  )dt  ). Since the pressure pin is an increasing function of time, to anticipate the

t

membrane permeability decrease the condition 0 pin (t  )dt  /t > 1 should hold, implying non-uniformity smaller than that of constant pressure case. Case 1b. Finite liquid recovery (finite value of A) In this case the pressure and velocity fields are described by the following system of ordinary differential equations (arising from the proper transformation of the second order equation for the pressure): dU −Ap = 1 + (1 + (p − 1))m dx

(38a)

dp = −U 1/(1+ε) dx

(38b)

The above system can be solved as initial value problem with U = Uin and p = pin (t) at x = 0. The deposit evolution is described by dm = uw dt p 1 + (1 + (p − 1))m

(39a)

The above equation must be solved numerically but in the particular case of  = 0 further analytical treatment is possible. In particular: Constant pressure

In the case of constant flux the function pin (t) must be computed to fulfill the requirement

m = −1 + (1 + 2(1 − x)t)1/2

1

uw = (1 + 2(1 − x)t)

−1/2

(35a)

uw =

uw dx = uwo,av

(35b)

The small t regime behavior (clean membrane) can be found through a first order expansion of the above relations, i.e. m = (1 − x)t, uw = 1 − (1 − x)t. Obviously, the fouling pattern follows the trans-membrane pressure pattern. In the large t limit, where the fouling contribution dominates the local membrane resistance, the following relations are obtained: m = (2(1 − x)t)1/2 , uw = (2(1 − x)t)−1/2 . The initial non-uniformity of deposit thickness is due to the nonuniform flux on the clean membrane; however, in regions where

(39b)

(40)

0

where uwo,av is the average permeate flux for the clean membrane (i.e. for m = 0). The variations of U and p along x are not very large so a simple Euler method is appropriate for solving at each time step. The time advancement is also made through the Euler method in the time domain. The numerical approach is fully explicit in the constant pressure case. However, these schemes cannot be directly used for the constant flux case where the condition in Eq. (40) must be

M. Kostoglou, A.J. Karabelas / Chemical Engineering Journal 187 (2012) 222–231

also satisfied. The problem can be overcome by adding the artificial equation

⎛

dpin = −E ⎝ dt

1

⎞

uw dx − uwo,av ⎠

(41)

equations, as suggested in the preceding Case (1b), results in the following system of equations for the complete problem: dU 1 = Ly dx

Ly

0

CASE 2. Finite pressure drop in the permeate side (finite value of B) Case 2a. Insignificant fluid recovery (A  1) As noted in Case 1a (ii), the pressure profile for constant inlet velocity can be described as p = pin − x, where  is constant in time. The governing equations are: ∂x2

+

uw =

∂2 P ∂y2

=

B(pin − x − P) 1 + (1 + (pin − x − P − 1))m

(pin − x − P) 1 + (1 + (pin − x − P − 1))m

(42) (43)

dm = uw dt

(44)

The numerical solution strategy in this case includes second order spatial discretization of the elliptic equation and solution of the resulting system at each time step through a modified Gauss–Seidel iterative scheme [31]. The P in the numerator and denominator of right hand side of Eq. (42) are treated in different ways (previous iteration values are used in denominator) in order to retain the linearity of the iterative scheme. The time advancement is done again through an Euler discretization of Eq. (44). At each time step, the Gauss–Seidel iterations for the pressure field start from the converged field of the previous step. This decreases the required number of iterations for convergence, and greatly accelerates the solution procedure. In case of constant flux strategy, the solution procedure can be retained by adding and discretizing properly the additional equation for pin inspired by the success of the corresponding treatment of the one dimensional problem (Eq. (41)):

⎛

dpin 1 = −E ⎝ Ly dt

−A(p − P) dy 1 + (1 + (p − P − 1))m

(46)

0

In this way and with the appropriate choice of parameter E, the explicit discretization of the mathematical problem can be retained even for the case of constant flux.

∂2 P

227

 

Ly 1

0

⎞

uw dxdy − uwo,av ⎠

(45)

0

Here uwo,av is the prescribed two-dimensional space average flux for the clean membrane. Case 2b. Finite flow recovery (finite value of A) This is actually the complete general case, based on the numerical approaches for the simpler cases and the results of formal mathematical analysis, for which a simplified numerical solution technique is developed. Kostoglou and Karabelas [8] showed that under practical conditions of membrane operation the complete pseudo-steady state problem can be expanded with respect to the parameter A to a series of linear problems. The implementation of this expansion procedure is not computationally advantageous in the case of two dimensional profiles of membrane permeability; however, this analysis is very important and helpful because it shows that the actual retentate profile is essentially one dimensional (with variation only in the x-direction) despite the two dimensional nature of membrane permeability and permeate pressure field. This result, therefore, guides further analysis. The governing equation for the one-dimensional retentate field is derived by taking the integral of Eq. (16) in the y direction and using the corresponding boundary conditions. Transformation of the resulting second order equation to a system of two first order

dp = −U 1/(1+ε) dx

(47)

∂2 P B(p − P) ∂2 P + = 1 + (1 + (p − P − 1))m ∂x2 ∂y2

(48)

uw =

(p − P) 1 + (1 + (p − P − 1))m

dm = uw dt

(49) (50)

The solution procedure is a combination of techniques suggested for the simplified cases including spatial Euler discretization of the initial value problem [Eqs. (46) and (47)], modified Gauss–Seidel iterative solution for the elliptic Eq. (48), and temporal Euler discretization of Eq. (50). In case of constant flux, solution of the additional Eq. (45) is needed. The overall solution procedure is by far more efficient compared to the direct numerical attack of the original problem. In order to get an idea of the difficulties of direct numerical attack to the system of Eqs. (16)–(28), it will be noted that even the pseudo-steady state problem with a fine finite element mesh constitutes a stiff problem and to achieve convergence takes significant time in a personal computer. To follow the temporal evolution of the deposit requires several decades of pseudo-steady state simulations. Further to this, in the case of constant flux an additional trial and error procedure must be employed at each time step to ensure constancy of total membrane flux. The alternative way of adding Eq. (45) to the system alleviates the need for the iterative procedure but increases the number of the required time steps leading to similar total computational time. 4. Indicative results and discussion Some indicative results regarding the distribution of the deposit mass m and of the flux uw on the membrane surface will be shown in the following. 4.1. No permeate side pressure drop (B = 0). At first, the case of no permeate side pressure drop with A = 0.15, Uin,av = 0.45 and ε = −0.4 is considered. These parameters correspond to initial values of R = 0.3 and po = 0.8 (i.e. 30% flow recovery, 20% pressure drop along the flow). The deposit mass and wall velocity profiles are always given normalized with their average values on the membrane in order to stress their non-uniformities. The evolution of the deposit mass and wall velocity profiles along the channel for the constant pressure case are shown in Figs. 2 and 3, respectively. The initial normalized profiles of deposit mass m and flux uw coincide and follow the clean membrane pressure profile, which is not linear in this case due to the finite A value. There are two important observations regarding the evolution of the m profiles. First, with increasing time (and correspondingly deposit growth) the profile tends to become linear. This is due to the fact that the average wall flux and correspondingly R tend to decrease due to fouling, so the pressure profile is getting linear as R → 0. Second, the reduction of the non-uniformity seems to be smaller (approx. 50%) than that derived in closed form in the preceding case for A  1 (see Eq. (35a)). The explanation for this is that, with time the wall flux tends to decrease and so the flow through the retentate channel increases leading to increased pressure drop (po

228

M. Kostoglou, A.J. Karabelas / Chemical Engineering Journal 187 (2012) 222–231

1.15

1.15

t=0 1.1

1.1 t=0.5 t=1

1.05

1.05

m/m av

av

m/m

t=2

1

t=5

1

0.95

t=0 t=1 t=5 t=50

0.95

0.9

0.9

0.85 0

0.2

0.4

0.6

0.8

0.85

1

x

0.2

0

Fig. 2. Evolution of normalized deposit mass profile (B = 0, A = 0.15, Uin,av = 0.45, ε = −0.4, constant pressure).

decreases in time). The increase of the pressure drop counterbalances the reduction of the non-uniformity from the stabilization effect of flux non-uniformity. The evolution of the wall velocity profiles in Fig. 3 reveals that they converge faster in time compared to the deposit profiles. The asymptotic profiles of m and uw are similar (as in the case of initial profiles) but this similarity does not hold for the intermediate profiles. The situation is quite different for the case of the constant mean flux membrane operation (i.e. R and Uin,av remain constant) for which the evolution of the m and uw normalized profiles are shown in Figs. 4 and 5, respectively; here the case of an incompressible cake (i.e.  = 0) is considered first. The initial profiles are by definition the same as in the case of constant pressure but their evolution is quite different. The increase of the inlet pressure decreases the fractional pressure drop leading to a fast reduction of non-uniformity. The decrease of non-uniformity is faster for the

0.4

0.6

0.8

1

x Fig. 4. Evolution of normalized deposit mass profile (B = 0, A = 0.15, Uin,av = 0.45, ε = −0.4, constant flux).

flux uw , which is expected since uw is an instantaneous quantity whereas m involves an integral in time thereby imparting an inertia feature in its dynamics. Figs. 4 and 5 provide an idea about the evolution of the non-uniformities for constant average wall flux, but they refer to a zero compressibility case; therefore, the effect of the compressibility must be also dealt with. Cake compressibility tends to enhance the non-uniformity reduction in two ways; first, it tends to reduce the local membrane permeability thereby increasing the regularizing effect of flux, and second it increases the inlet pressure, thus reducing the fractional pressure drop. The increase of the inlet pressure for the specific example presented here, for several values of the compressibility parameters, is depicted in Fig. 6. It is interesting that in the cases depicted in Figs. 2–5 the profiles evolve in time by maintaining a constant value at a fixed point in the axial direction. 1.15

1.15

t=0 1.1

1.1

t=0.2

t=0 t=1 t=5 t=10

t=0.5 1.05 t=1

u /u

w

1

w

u /u

w,av

w,av

1.05

t=2 1

0.95

0.95

0.9

0.9

0.85 0

0.2

0.4

0.6

0.8

1

x Fig. 3. Evolution of normalized permeate flux profile (B = 0, A = 0.15, Uin,av = 0.45, ε = −0.4, constant pressure).

t=5

0.85 0

0.2

0.4

x

0.6

0.8

1

Fig. 5. Evolution of normalized permeate flux profile (B = 0, A = 0.15, Uin,av = 0.45, ε = −0.4, constant flux).

M. Kostoglou, A.J. Karabelas / Chemical Engineering Journal 187 (2012) 222–231

3

229

1.3 t=0

x=0

=0 =0.1 =0.2 =0.3

2.5

t=0.2 1.2 t=0.5 t=1 t=2

w,av

1.1

t=3

u /u

pin

t=4

w

2

1

x=1 t=3

1.5

0.9

t=4 t=2 t=1

t=0.5 t=0.2

t=0 0.8 0

1 0

0.5

1

1.5

0.2

0.4

2

t

0.6

y

0.8

1

Fig. 8. Evolution of normalized permeate flux profile (B = 0.5, A = 0, po = 0.8, constant pressure) in the lateral (y) direction for x = 0 and x = 1.

Fig. 6. Effect of fouling layer compressibility parameter  on inlet pressure evolution (B = 0, A = 0.15, Uin,av = 0.45, ε = −0.4, constant flux).

4.2. Pressure drop at permeate side considered To demonstrate the influence of the permeate side pressure drop on the non-uniformity of the deposit, the case of A = 0, B = 0.5, po = 0.8 for constant pressure operation is considered. The evolution of the normalized m and uw profiles (in the y direction at the inlet (x = 0) and outlet (x = 1) of the retentate channel) is presented in Figs. 7 and 8, respectively. The curves do not have an average value of unity because the normalization refers to the complete 2-dimensional area of the membrane. The deposit mass is higher at the inlet of the retentate channel due to the greater value of the pressure there. Also, the mass is greater close to the permeate flow outlet (y = 0) due to the lower permeateside pressure which leads to an increase of the trans-membrane

pressure there. It is interesting that at x = 0 mainly the small y part of the m profile evolves, whereas at x = 1 mainly the large y part evolves. Obviously the degree of non-uniformity decreases in time in both x and y directions. Similar profiles are shown for the permeate flux in Fig. 8, the only difference being their larger changes. The evolution of the permeate side pressure along the line y = 1 where the permeate pressure is maximum is shown in Fig. 9. With passing time the wall flux decreases due to fouling, so the pressure drop and correspondingly the absolute pressure in the permeate side is reduced. In the particular case examined here, the maximum pressure in the permeate side is at approx. 18% of the operating pressure, as shown in Fig. 9, which can be considered large for practical membrane operations. This value is used here to better stress the non-uniformities induced by transmembrane pressure differences. The permeate pressure displays

0.2 1.3 t=0

t=0

x=0

0.18

t=0.2 1.2

t=0.5 t=1

t=0.2

0.16

t=0.5

P 0.14

m/m

av

t=2 1.1

t=1 0.12

1

x=1 t=2

0.1 0.9

t=0

0.8 0

0.2

0.4

y

0.6

0.8

t=3

t=1 t=0.5 t=0.2

0.08

1

Fig. 7. Evolution of normalized deposit mass profile (B = 0.5, A = 0, po = 0.8, constant pressure) in the lateral (y) direction for x = 0 and x = 1.

0

0.2

0.4

0.6

0.8

1

x Fig. 9. Evolution of permeate side pressure along y = 1 line (B = 0.5, A = 0, po = 0.8, constant pressure).

230

M. Kostoglou, A.J. Karabelas / Chemical Engineering Journal 187 (2012) 222–231

(a)

1.3 1.2 1.1 1 0.9 M 0.8

0 0.2 0.4 0.8

x

1 0.2

0

0.4

0.6

0.6

3(1 − /2)

y

(b)

1.3 1.2 1.1 1 0.9 M 0.8 0.7

0 0.2 0.4 0.6

x

0.6 0.5 0.8 1 0

0.2

0.4

0.6

0.8

1

y

Fig. 10. Two dimensional profile of normalized deposit mass M = m/mav for B = 0.5, A = 0, po = 0.8 at (a) t = 0 and (b) t = 2.

a small variation along the x-direction, which tends to vanish with time. Although the 1-dimensional profiles at several (x or y) positions on the membrane are informative to reveal important features of the evolution of the deposit, 2-dimensional ones are needed to provide a complete picture of their form at a certain point in time. Such profiles for clean and fouled membranes (A = 0, B = 0.5, po = 0.8) are shown in Fig. 10a and b, respectively. These graphs clearly show the previously discussed initial non-uniformity of the mass density distribution throughout a membrane sheet, and its subsequent reduction (“smoothing”) with time, which appears to be a notable outcome of the dynamic interaction of the deposition process with the flow fields at both sides of the membrane. From the above results it is clearly evident that the pressure profile tends to create (modestly) non-uniform fouling on the membrane. A key question is how the cake non-uniformity influences the integral properties of the membrane, which are of practical interest, i.e. the evolution of the fluid recovery R for the case of constant inlet flow rate and pressure, and the evolution of pressure for constant inlet flow rate and mean flux R. A first answer to these questions can be obtained by reformulating the mathematical model, assuming a uniform distribution of the cake throughout the membrane; thus,

(51)

resulting from integration of Eq. (18) on the membrane surface. The comparison of problem solutions with spatially distributed m and with uniform m reveals a relative error for R and pin of the order of 1/1000. This small error can be shown more clearly by considering the case treated by Eq. (34). The ratio of approximate to exact average mass for t  1 can be shown to be 2(1 − (1 − )3/2

1

0.8

dmav = uw,av dt

0.7

0.5

0.6

a spatial average mav replaces the function m(x,y) in the governing equations. The evolution of mav is described by the equation

1/2

.

This expression becomes smaller than 0.999 for  > 0.26 which means that it is totally negligible for any practical value of pressure drop. Therefore, the cake non-uniformity does not seem to have an effect on the integral membrane properties. It should be noted that the uniform cake approximation does not directly simplify at all the mathematical problem that must be solved (since the structure of the problem remains the same); however, such approximate solutions for the case of uniform membrane permeability have already been developed for computational convenience [5,32]. Indeed, these approximate solutions can be combined with the cake evolution equation, to tackle the complete problem with no need to solve differential equations in the spatial domain. In conclusion, the organic–gel cake developing on flat-sheet membranes (similar to those comprising SWM), in the case of weak ionic strength feedwaters (effectively neglecting the presence of salts), exhibits spatial non-uniformities. In the case of non-zero pressure drop in the permeate side, these non-uniformities exist in both axial and lateral (x and y) directions. Only the x-direction non-uniformities develop for the case of negligible pressure drop at the permeate side. The x-direction non-uniformities depend on the pressure drop along the retentate channel. The results obtained by considering the dynamic interactions of species deposition with the flow fields, show that the deposit formation process tends to be “self-stabilized” (leading to a kind of deposit smoothing) in that the initial (normalized) non-uniformities arising from the pressure non-uniformities tend to be reduced as the cake grows, but they are not completely eliminated. The reduction of non-uniformities appears to be faster for the constant flux than for the constant pressure membrane operation. By increasing the fluid recovery fraction R, the rate of reduction of the non-uniformities tends to decrease. Finally, the fouling layer non-uniformities themselves (under the conditions studied) do not appear to influence the overall membrane operation features, which are of course affected by mean deposit mass density and the related added resistance to fluid permeation. It is noted here that this result is obtained for the relatively simple (essentially linear) deposit growth mode for the particular organic/colloidal deposits, but it cannot be generalized for any kind of deposit. For example, for some cases of solid deposits (scale) created due to concentration polarization of sparingly soluble salts (like calcium carbonate), the deposition mechanism is extremely non-linear and even small differences of permeate flux can lead to large differences of deposition rate [33]; the results of the present work are not applicable to such cases. 5. Conclusions A theoretical study is presented regarding the spatial and temporal evolution of deposits and permeate flux on a membrane

M. Kostoglou, A.J. Karabelas / Chemical Engineering Journal 187 (2012) 222–231

sheet (similar to that comprising spiral wound modules) for the case of a linear deposit growth mechanism, observed to be typical of organic fouling layers. Fouling layer (cake) compressibility is introduced through a rather simple constitutive relation to facilitate mathematical analysis. The mathematical model is formulated in terms of specific parameters that determine the dimensionality of the solution. Several analytical and asymptotic results are derived for the one-dimensional case expanding the zerodimensional results that are in agreement with well known results for dead-end filtration [1]. The two cases of constant inlet pressure and constant average flux (representing practically realistic modes of operation) are treated separately. Unlike the zerodimensional case, the constant total flux mode does not imply constant local flux since only the average flux is kept constant, with the local flux distribution evolving in time. Efficient numerical methods have been developed hierarchically from one to two dimensions and they can serve as the basis for tackling more complicated deposit formation mechanisms in the future. The results, obtained here for the specific linear deposit growth mechanism, suggest that ignoring the deposit spatial non-uniformity has an insignificant effect on the computation of the macroscopic process parameters. References [1] Water Treatment Principles and Design, Chapter 12, Membrane Filtration in Water Treatment, Wiley, New York, 2005. [2] J. Schwinge, P.R. Neal, D.E. Wiley, D.F. Fletcher, A.G. Fane, Spiral wound modules and spacers. Review and analysis, J. Membr. Sci. 242 (2004) 129–153. [3] C.P. Koutsou, S.G. Yiantsios, A.J. Karabelas, Direct numerical simulation of flow in spacer-filled channels: effect of spacer geometrical characteristics, J. Membr. Sci. 291 (2007) 53–69. [4] C.P. Koutsou, S.G. Yiantsios, A.J. Karabelas, A numerical and experimental study of mass transfer in spacer-filled channels: effects of spacer geometrical characteristics and Schmidt number, J. Membr. Sci. 326 (2009) 234–251. [5] M. Kostoglou, A.J. Karabelas, On the fluid mechanics of spiral-wound membrane modules, Ind. Eng. Chem. Res. 48 (2009) 10025–10036. [6] D. Sioutopoulos, S.G. Yiantsios, A.J. Karabelas, Relation between fouling characteristics of RO and UF membranes in experiments with colloidal organic and inorganic species, J. Membr. Sci. 350 (2010) 62–82. [7] D. Sioutopoulos, A.J. Karabelas, S.G. Yiantsios, Organic fouling of RO membranes: Investigating the correlation of RO and UF fouling resistances for predictive purposes, Desalination 261 (2010) 272–283. [8] M. Kostoglou, A.J. Karabelas, Mathematical analysis of the meso-scale flow field in spiral wound membrane modules, Ind. Eng. Chem. Res. 50 (2011) 4653–4666. [9] D. Hasson, A. Drak, R. Semiat, Inception of CaSO4 scaling on RO membranes at various water recovery levels, Desalination 139 (2001) 73–81. [10] C. Tzotzi, T. Pahiadaki, S.G. Yiantsios, A.J. Karabelas, N. Andritsos, A study of CaCO3 scale formation and inhibition in RO and NF membranes, J. Membr. Sci. 296 (2007) 171–184.

231

[11] A.I. Schafer, A.G. Fane, T.D. Waite, Nanofiltration of natural organic matter: removal, fouling and the influence of multivalent ions, Desalination 118 (1998) 109–122. [12] S. Lee, M. Elimelech, Relating organic fouling of reverse osmosis membranes to intermolecular adhesion forces, Environ. Sci. Technol. 40 (2006) 980–987. [13] Y. Ye, P. Le Clech, V. Chen, A.G. Fane, B. Jefferson, Fouling mechanisms of alginate solutions as model extracellular polymeric substances, Desalination 175 (2005) 7–20. [14] J.S. Vrouwenvelder, D.A. Graf von der Schulenburg, J.C. Kruithof, M.L. Johns, M.C.M. van Loosdrecht, Biofouling of spiral-wound nanofiltration and reverse osmosis membranes: a feed spacer problem, Water Res. 43 (2009) 583–594. [15] J.C. Chen, A.S. Kim, Monte Carlo simulation of colloidal membrane filtration: principal issues for modeling, Adv. Colloid Interface Sci. 119 (2006) 35–53. [16] S.G. Yiantsios, D. Sioutopoulos, A.J. Karabelas, Colloidal fouling of RO membranes: an overview of key issues and efforts to develop improved prediction techniques, Desalination 183 (2005) 257–272. [17] F. Evangelista, Optimal design and performance of spiral wound modules. I. Numerical method, Chem. Eng. Commun. 72 (1988) 69–81. [18] V.V. Ranade, Computational Flow Modeling for Chemical Reactor Engineering, Academic Press, New York, 2002. [19] Y.S. Polyakov, Deadend outside-in hollow fiber membrane filter: mathematical model, J. Membr. Sci. 279 (2006) 615–624. [20] R. Aris, Mathematical Modeling Techniques, Dover, New York, 1994. [21] V.V. Ranade, A. Kumar, Comparison of flow structures in spacer filled flat and annular channels, Desalination 191 (2006) 236–244. [22] G.A. Fimbres-Weihs, D.E. Wiley, Review of 3D CFD modeling of flow and mass transfer in narrow spacer field channels in membrane modules, Chem. Eng. Process. 49 (2010) 759–781. [23] S. Wardeh, H.P. Morvan, CFD simulations of flow and concentration polarization in spacer filler channels for application to water desalination, J. Chem. Eng. Res. Des. 86 (2008) 1107–1116. [24] C.P. Koutsou, S.G. Yiantsios, A.J. Karabelas, Numerical simulation of the flow in a plane channel containing a periodic array of cylindrical turbulence promoters, J. Membr. Sci. 231 (2004) 81–90. [25] G. Guillen, E.M.V. Hoek, Modeling the impacts of feed spacer geometry on reverse osmosis and nanofiltration processes, Chem. Eng. J. 149 (2009) 221–231. [26] E.M.V. Hoek, M. Elimelech, Cake-enhanced concentration polarization: a new fouling mechanism for salt rejecting membranes, Environ. Sci. Technol. 37 (2003) 5581–5588. [27] B.F. Ruth, Correlating filtration theory with industrial practice, Ind. Eng. Chem. 328 (1946) 564–571. [28] E. Vorobiev, Derivation of filtration equations incorporating the effect of pressure redistribution on the cake-medium interface: a constant pressure filtration, Chem. Eng. Sci. 61 (2006) 3686–3697. [29] C. Tien, R. Bai, An assessment of the conventional cake filtration theory, Chem. Eng. Sci. 58 (2003) 1323–1336. [30] S. Sundaramoorthy, G. Srinivasan, D.V.R. Murthy, An analytical model for spiral wound reverse osmosis membranes modules: Part I – model development and parameter estimation, Desalination 280 (2011) 403–411. [31] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, New York, 1992. [32] F. Evangelista, An improved analytical method for the design of spiral wounds modules, Chem. Eng. J. 38 (1988) 33–40. [33] A.J. Karabelas, M. Kostoglou, S.T. Mitrouli, Incipient crystallization of sparingly soluble salts on membrane surfaces: the case of dead-end filtration with no agitation, Desalination 273 (2011) 105–117.