- Email: [email protected]
Simulation of fouling and backwash dynamics in dead-end microfiltration: Effect of pore size
Journal of Membrane Science 392–393 (2012) 48–57
Contents lists available at SciVerse ScienceDirect
Journal of Membrane Science journal homepage: www.elsevier.com/locate/memsci
Simulation of fouling and backwash dynamics in dead-end microfiltration: Effect of pore size Tsutomu Ando a,∗ , Kazuki Akamatsu b , Shin-ichi Nakao b , Masahiro Fujita c a b c
Department of Mechanical Engineering, College of Industrial Technology, Nihon University, 1-2-1 Izumi-cho, Narashino, Chiba 275-8575, Japan Department of Environmental and Energy Chemistry, Kogakuin University, 2665-1 Nakano-machi, Hachioji, Tokyo 192-0015, Japan Department of Chemical System Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
a r t i c l e
i n f o
Article history: Received 2 September 2011 Received in revised form 29 November 2011 Accepted 30 November 2011 Available online 9 December 2011 Keywords: Membrane pore size Particle fouling Dead-end microfiltration Backwash Simulation
a b s t r a c t We carried out numerical simulations of pressure-driven dead-end microfiltration with a backwash operation using a two-way coupling model with consideration of the particle–fluid interactions. This study numerically modeled a membrane with regularly spaced straight pores, which are assumed to be tracketched pores. On the basis of the results obtained by microfiltration under fixed particle concentration (5%), the present paper examines the effect of pore size on particle fouling by comparing snapshots of particle motion, the permeate flux, and the total resistance of the membrane over a period of time. The results of simulations confirm two modes for particle fouling when particle diameter d = 100 nm. The larger pore size membrane (dm = 3.6d) shows a fouling mode in which initially particles filled the pores and then form a cake layer on the surface of membrane. The smaller pore size membrane (dm = 2.5d) shows a fouling mode in which particles are accumulated on the surface of the membrane without filling the pores and a cake layer forms across the entire filtration. These general behaviors are in agreement with the hypothesis deduced from experiments. In addition, we carried out simulations of a sequence of filtration including a backwash operation and compared the fouling conditions for filtration before and after its inclusion. The particles remaining on the membrane after the backwash operation decrease the effective pore area and cause a change in fouling conditions from pore blockage to cake formation. © 2011 Elsevier B.V. All rights reserved.
1. Introduction In recent years, attention has been given to water treatment using membrane filtration, e.g., water purification, wastewater treatment and seawater desalination [1–3]. In general water treatment, nanofiltration and reverse osmosis primarily target the separation of ions, ultrafiltration targets the separation of molecules, and microfiltration targets the separation of particles and microorganisms. One of the most significant contemporary issues in microfiltration in the water treatment process is a decrease in membrane performance through fouling caused by particles suspended in water. As such, estimation of the decrement in microfiltration is essential in optimizing the membrane filtration process [4,5]. However, before arguing this issue, we should first and foremost understand the mechanism of fouling in microfiltration. What is the origin of fouling, and what are the precursors to said origin? There are many things we do not yet understand regarding fouling in microfiltration.
∗ Corresponding author. Tel.: +81 47 474 2338; fax: +81 47 474 2349. E-mail address: [email protected] (T. Ando). 0376-7388/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.memsci.2011.11.051
Since ultrafiltration and reverse osmosis are controlled by the diffusive motions of molecules, their performance can be partially estimated using the molecular diffusive model [6]. On the other hand, the hydrodynamic effect is dominant in microfiltration because particle motion is affected by solvent flow. Therefore, simulation of microfiltration requires the fluid dynamic model rather than the molecular diffusive model. Recently, many computational fluid dynamics (CFD) simulations have been used to evaluate the microfiltration process [7–9]. However, most of these extant CFD filtration simulations were based on the one-way method. This method first resolves the flow field around the membrane, and then particles are traced along the streamline in the flow field. In this case, particle motion does not affect solvent flow around the particles themselves; therefore, the one-way method cannot precisely resolve questions of concentrated particulate flow. The authors have worked toward developing a meso-scale model in which we consider the dynamics of Brownian particles and the electrostatic force between two particles. We have also conducted simulation studies of meso-scale concentrated particulate flow. This simulator, called “SNAP-F (Structure of NAnoParticles Flow),” is a two-way coupling model that can simultaneously resolve both the motions of particles and fluid while also incorporating particle–fluid interactions [10–12]. As the lattice spacing
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
49
used to resolve the flow field is much smaller than the diameter of particles and the flow of the surface on particles is resolved with high resolution, we can calculate the solvent flow through the fouling particles. SNAP-F applied to the microfiltration process can thereby estimate the filtration performance of a membrane system. Moreover, in addition to the process of surface deposition and the fouling of pores, SNAP-F can also simulate the backwash dynamics for the detachment of particles from both the surface of the membrane and the inside of the pores; this reveals the quantitative effect of the actual filtration process including both fouling and backwashing. Previously, in two-dimensional simulations of concentrated colloidal nanoparticles, we have shown that our model of solid contact interactions, including frictional forces between nanoparticles, describes the dynamics of aggregated colloidal nanoparticles [10], and we have discussed the aggregation and dispersion structure of nanoparticles (depending on the Peclet number) in simple shear flows [11]. In a previous paper [12], we developed a microfiltration membrane model with a single straight pore and carried out three-dimensional simulations of concentrated particulate flow in pressure-driven dead-end microfiltration. We discussed the filtration performance to show the permeate flux and particle rejection with visualization of particle fouling on a membrane model with a single straight pore. The present paper carries out simulations of pressure-driven dead-end microfiltration with a backwash operation. Here, we apply the membrane model with four regularly spaced cylindrical pores, which are assumed to be track-etched pores. This study addresses three kinds of membranes with more than twice and less than four-times pore size of the particle diameter. Since these pore diameters are larger than the particle diameter and are not large enough to accumulate on the inside wall of the pore, fouling is not classified by the conventional fouling model (such as “complete blocking” and “standard blocking”) [13,14]. On the basis of the results of microfiltration, we examine the effect of pore size on particle fouling by comparing snapshots of particle fouling, the permeate flux and the total resistance of membrane over a period of time. We then perform backwash simulation for the membrane with pore size three-times that of the particle diameter. Based on the results of the backwash operation, we discuss the fouling conditions for filtration before and after the backwash operation.
In the present paper, we used a two-step method [17,18]. The previous paper used a one-step method and fixed the pressure in an inner region of particles to zero in order to satisfy the conservation of mass on fluid–particle coupling, but the equation of continuity was not satisfied on the boundary of the particles. This study divides the pressure into two to satisfy the conservation of mass: the intermediate pressure without particles, p∗ , and the correcting pressure with particles, p .
2. Simulation model
∇ · vn+1 = 0
2.1. Flow of solvent
Then, operating on both sides of Eq. (6), and substituting it for Eq. (8), the Poisson equation of intermediate pressure is given by
The momentum equation of the entire flow field, containing both particles and solvent is the fluctuating Navier–Stokes equation [15], shown as follows [12].
∇ 2 p∗ =
∂v 1 1 1 + (v · ∇ )v = − ∇ p + ∇ 2 v + ∇ · S − D + ␣ ∂t ␣=
vp − v 1 1 + (v · ∇ )v − ∇ 2 v − ∇ · S + D t
p = p∗ + p
(3)
Substituting Eq. (3) for Eq. (1), we obtain the following equation using the flow velocity without particles v∗ .
1 1 1 ∂v∗ = (1 − ) −(v∗ · ∇ )v∗ + ∇ 2 v∗ + ∇ · S − D − ∇ p∗ ∂t +
vp − v∗ t
−
1 ∇ p
(4)
Here the acceleration of solvent received from the rigid body is defined by ␣ =
vp − v∗ t
+ (v∗ · ∇ )v∗ −
2 ∗ 1 1 1 ∇ v − ∇ · S + D − ∇ p (5)
Eq. (4) is divided into 2 steps predictor–corrector method as follows. (1st step)
v ∗ − vn t
= (1 − ) −(v∗ · ∇ )v∗ +
(2nd step)
vn+1 − v∗ = t
according
to
the
1 2 ∗ 1 1 ∇ v + ∇ · S − D − ∇ p∗
vp − v∗ t
−
1 ∇ p
(6)
(7)
For the 1st step, the continuity equation without particles is expressed by
∇ · v∗ = 0
(8)
For the 2nd step, the continuity equation of the new entire flow field with particles is expressed by
(9)
∇ · vn t
+∇ ·
(1 − )
−(v∗ · ∇ )v∗ +
2 ∗ 1 1 ∇ v + ∇ ·S− D
(10)
(1)
Then also, operating on both sides of Eq. (7), and substituting it for Eq. (9), the Poisson equation of correcting pressure is given by
(2)
∇ 2 p =
where t denotes time, v is the fluid velocity vector on the computational lattice, is the solvent density, p is the pressure, is the solvent viscosity, S is the thermal fluctuating stress tensor, and D is the pressure gradient vector. is the volume fraction of particles and membrane, and this value is [0,1] considering the overlap between two particles or between a particle and the membrane. ˛ is the acceleration vector associated with the rigid velocity of particles on the lattice. Our previous paper [12] defined the membrane using the FSA (fortified solution algorithm) [16], but the present paper treats the membrane the same as particles, i.e., the FSA is not used. Here the velocity of the membrane is fortified as zero.
∇ · v∗ t
vp − v∗
+∇ ·
t
(11)
The algorithm to evolve the entire flow field from n-step to n + 1step is given as follows: 1. Solve the prediction values v∗ , p∗ from the present value vn using Eqs. (6) and (10). 2. Solve the new value vn+1 and the correction value p from the prediction value v∗ using Eqs. (7) and (11). 3. Calculate the new value pn+1 from the prediction value p∗ and the correction value p using Eq. (3). 4. Calculate the force and torque acting on particles using ␣ of Eq. (5).
50
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57 Table 1 Simulation condition.
2.2. Motion of particles The motion of particles in solvent is governed by Newton’s equation of motion. The translational motion of the l-th particle is expressed by the following equation [12]. ml
dV l = F cl + F el + F vl + F hl + F lb l dt
Particle
(12)
where ml is the mass of the particle, Vl is the translational velocity vector of the particle, Fcl is the contact force, Fel is the electrostatic v
force, Fl is the van der Waals force,
Fhl
is the hydrodynamic force,
and Flb l is the lubrication force. The rotational motion of the l-th particle is expressed by the following equation. Il
dl = Tcl + Thl dt
(13)
where Il is the moment of inertia of the particle, l is the angular velocity vector of the particle, Tcl is the contact torque, and Thl is the hydrodynamic torque. The contact force and torque are calculated using the Voigt model and by considering the Hertzian theory and the Mindlin model [10,11,19–22]. Here, the normal and the tangential contact forces both between two particles and between a particle and the membrane are given by a combination of elastic force and damping force. Then, the electrostatic force and the van der Waals force both between two particles and between a particle and the membrane are given by the Derjaguin–Landau–Verwey–Overbeek (DLVO) theory [23]. In addition, as in the previous paper [12], the present study considers the normal component of the lubrication force between pairs of particles or between a particle and a plane wall when the gap width tends to zero [24,25]. Because the lubrication force is a dominant force on a scale that is smaller than the lattice resolution, the analytical expressions of the lubrication force both between particles and between a particle and the membrane have been added to compensate the underestimation by the lack of lattice resolution in this study. The hydrodynamic force F hl and torque Thl acting on the l-th particle are given as below, using the volume integration of ␣ .
F hl
l (r){␣ (r) + D}dV
(14)
l (r){r l (r) × ␣ (r)}dV
(15)
p
=− Vp
T hl = −
p
Vp
Solvent
where V p denotes the region occupied by the l-th particle. The diffusion of particle is obtained through the hydrodynamic force F h and torque T h including the stochastic fluctuating stress. This results in the many-body Brownian motion of the particles.
Membrane
Nomenclature
Value
Diameter,d [nm] Density,p [kg/m3 ] Young’s modulus,E p [Pa] Poisson ratio, [–] Coefficient of friction between particles, p [–] Coefficient of friction between particle and membrane, m [–] Zeta potential, p [mV] Concentration, v [%] Hamaker constant, A [J]
100 3.97 × 103 3.65 × 109 0.254 0.1
Density, [kg/m3 ] Viscosity, [Pa s] Relative permittivity, εr [–] Ionic strength of the electrolyte, I [mol/l] Valence of ions [–] Temperature, T [K]
1.0 × 103 1.0 × 10−3 80 10−4
Thickness, t m [nm] Pore diameter, dm [nm] Porosity, ε [–] Density, m [kg/m3 ] Young’s modulus, Em [Pa] Poisson ratio, m [–] Zeta potential, m [mV]
400 250, 300, 360 0.306, 0.441, 0.636 3.97 × 103 3.65 × 109 0.254 −40
0.1 −40 5 4.7 × 10−20
1 293.15
filtration region, as shown in Fig. 1(a), because we employ a periodic boundary condition on the faces parallel to the direction of flow. The computational domain is (x, y, z) = (16d, 8d, 8d), as shown in Fig. 1(b). The membrane is 4d (400 nm) in thickness with 4 regularly spaced straight cylindrical pores of standard diameter 3d (300 nm) and is located in the middle of the flow direction of the computational domain. This study has simulated the real microfiltration, but the membrane is thinner than the real size due to limitations of simulation time. The lattice spacing of the flow field is 1/9d. As the initial condition, there are some particles in the feed-side domain according to the set value of particle concentration. After the simulation begins, the flow is induced by a uniform pressure gradient, and particles are generated at a point in the y–z plane on the feed side using random numbers, depending on the particle concentration. At the boundaries on the feed and permeate side, which are perpendicular to the direction of flow, the Neumann condition ıvx /ıx = 0 is employed and vy = vz = 0. As for the time step, we set the value as t = 0.4 ns so as not to induce numerical instability in all cases. 3. Simulation results and discussion
2.3. Simulation method and conditions
3.1. Comparative simulation between 1 pore and 4 pores
The equation of translational and rotational motions of particles is solved by using the first-order Euler explicit scheme, and the equation of particle position is solved by using the Crank–Nicolson scheme. Then, the fluctuating Navier–Stokes equation is solved using the two-step method; the first step uses the Semi-implicit method for pressure-linked equations shortened (SIMPLEST) [26], and the second step uses the convergence algorithm by the iterative calculation of pressure. Table 1 shows the condition of filtration simulation in the present study, and Fig. 1 shows the computational domain. The solvent is assumed to be water, and the particle is a spherical body with a diameter d (d = 100 nm) with physical properties identical to those of an alumina particle. In the simulation of pressure-driven dead-end microfiltration, the computational domain is a part of the
Without particles (v = 0%), we carried out simulations of pressure-driven dead-end microfiltration of pore size dm = 3d (300 nm), which is an ideal track-etched pore. The pressure gradient Dx = −2 × 1011 Pa/m is applied due to limitation of simulation time, and the difference in pressure between the feed and the permeate surfaces of the membrane is about 320 kPa. Then, the solvent flow becomes the Hagen–Poiseuille flow in the pore and the maximum Reynolds numbers Rep (= d|vx |)/ at both S = 0 and S = / 0 are around 0.42 for both the 1-pore membrane and the 4-pore membrane under the same porosity (ε = 0.441). The flow accelerates just before the tracer particles reach the entrance of the pore and flow inward, decelerates rapidly at the pore exit, and diffuses radially in the permeate domain, as seen in the figure without particles in Ref. [12], where simulation of microfiltration
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
51
Fig. 2. Permeate flux as a function of time up to t = 10 s for the 1-pore membrane and the 4-pore membrane, pore diameter dm = 3d.
flux is given by the average flow rate in the pore divided by the cross-sectional area of membrane (8d*8d). As the permeate volume increases, the fluxes of all pore sizes become smaller and indistinguishable due to resistance from particle fouling becoming greater than the resistance of the membrane. The permeate flux data can be used to examine the fouling condition by replotting the data in terms of the total resistance Rt , as given by the following equation: Rt =
Fig. 1. Schematic of calculation domain for studying dead-end microfiltration: (a) periodic boundary condition on the parallel faces with respect to the direction of flow and (b) computational domain, including membrane, with four regularly spaced straight pores.
was performed using the membrane model with a single straight pore. Fig. 2 shows the permeate fluxes, J, of each pore obtained by two simulations of the 1-pore membrane and the 4-pore membrane of pore size dm = 3d (ε = 0.441), under a particle concentration of v = 5%. Here, the permeate flux J of each pore for the 4 pore is given by the average flow rate in the pore divided by the 1/4 cross-sectional area of membrane (4d*4d). As each pore for the 4pore membrane draws an independent flux curve, simulation with as many pores as possible is important to obtain the image of the whole membrane, as shown in Fig. 1(a). Therefore, this study uses the membrane model with 4 regularly spaced straight pores and discusses the effect of pore size on particle fouling in microfiltration.
P J
(16)
The results are shown in Fig. 4, where P is the transmembrane pressure (320 kPa). The total resistance for the smaller pore size membrane of dm = 2.5d is initially larger due to the smaller pores but only increases with a decreasing slope for the entire simulation. The total resistance of the larger pore size membrane of dm = 3d and 3.6d starts off lower but becomes larger than that of the smaller pore size membrane by the end of the microfiltration. The total resistance curves of the larger pore size membrane have a point of inflection and produce an S-shape. These general behaviors agree with the experimental results of Tracy and Davis [27].
3.2. Simulation of filtration for different pore sizes We carried out simulations of dead-end microfiltration for two kinds of pore size, dm = 3.6d (360 nm) and 2.5d (250 nm), under a fixed particle concentration of v = 5%. The porosities of each pore size membrane are ε = 0.636 and 0.306, and 1.5 and 1/1.5 times that of standard pore (dm = 3d, ε = 0.441), respectively. Regarding the three kinds of pore sizes, we compared the permeate flux versus total permeate volume, as shown in Fig. 3. Here the permeate
Fig. 3. Permeate flux as a function of permeate volume for three kinds of pore diameter dm .
52
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
Fig. 4. Total resistance as a function of time for three kinds of pore diameter dm .
We first observed snapshots of the fouling condition of the standard pore size dm = 3d, as shown in Fig. 5. The total resistance curve for dm = 3d increases with an increasing slope to the point of inflection (at around 5 s), after which resistance continues to increase but with a decreasing slope. Compared with Fig. 5(a)–(d), which are snapshots at 4, 5, 6, and 7 s, respectively, it can be seen that pores plug at 6 s. We can also confirm this fact from Fig. 6, which shows the corresponding probability of particle existence. t = 6 s is the time when the flux curve for dm = 3d intersects with that for dm = 2.5d in Fig. 3, that is, the flux J is 0.16 m3 /(m2 s). However, the fouling mode does not strictly alternate to either pore blockage or cake formation at t = 6 s, and, in fact, it is assumed that the cake formation starts from around t = 5 s of the inflection point. Next, we examined two simulations of dm = 3.6d as a lager pore and 2.5d as a smaller pore to compare the total resistance versus time curve. The snapshots of dm = 3.6d and 2.5d are shown in Figs. 7 and 8, respectively. For the larger pore size membrane (dm = 3.6d), the positive second derivative of the total resistance curve initially fouls pores internally, and then the negative second derivative of the total resistance curve forms the cake layer on the surface of the membrane, as shown in Figs. 4 and 7, respectively. This behavior is consistent with the discussions of Tracy and Davis [27] and Ho and Zydney [28]. For the smaller pore size membrane (dm = 2.5d), the second derivative of the total resistance curve is negative for the entire simulation, and the fouling appears to be due to a growing cake layer on the surface of the membrane, as shown in Figs. 4 and 8, respectively. This behavior is also consistent with the discussions of Tracy and Davis [27]. In addition, we quantitatively compare the time transition of the fouling condition to examine the corresponding probability of particle existence of dm = 3.6d and 2.5d, as shown in Figs. 9 and 10, respectively. As shown in Fig. 9, in the case of dm = 3.6d, by 1 s, the probability of particle existence between the entrance and the exit of the pore is already at nearly the same value, 0.22 and 0.23. On the other hand, at this time, the probability of particle existence inside the pore is less than half that of the entrance and the exit. At t = 2 s, the particle existence at the pore exit is a 79% increase from 1 s, but at the same time at the entrance, there is only a 15% increase, and that at the inside center of the pore (x = 8d) has changed very little. The probability of particle existence at t = 5 s draws a monotonically decreasing curve from the exit to the entrance of the pore. After t = 10 s, the probability of particle existence only increases in the feed side of the membrane. In summary, for ideal, numerically modeled straight cylindrical pores, the fouling of dm = 3.6d happens initially at the exit and thereafter
Fig. 5. Snapshots of simulation for pore diameter dm = 3d: (a) 4 s, (b) 5 s, (c) 6 s, and (d) 7 s.
propagates to the entrance and inside of the pore, as shown in Fig. 7. In other words, in the early stages of filtration without particles in the pore, the Poiseuille flow in the pore decelerates rapidly at the exit and diffuses radially in the permeate domain. Then, in tandem with this solvent flow, some rolling particles on the inside wall remain on the permeate-side surface of the membrane when
Fig. 6. Probability of particle existence for pore diameter dm = 3d.
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
53
Fig. 7. Snapshots of simulation for pore diameter dm = 3.6d: (a) 1 s, (b) 2 s, (c) 5 s, (d) 10 s, and (e) 15 s.
Fig. 8. Snapshots of simulation for pore diameter dm = 2.5d: (a) 1 s, (b) 2 s, (c) 5 s, (d) 10 s, and (e) 15 s.
particles vacate the pore. These remaining particles around the exit reduce the effective cross-sectional area of the pore and become the origin of fouling. As shown in Fig. 10, in the case of dm = 2.5d, at t = 1 s, the probability of particle existence at the pore exit is less than 0.1 and increases to 0.29 until t = 5 s; thereafter, it does not increase. Conversely, the probability of particle existence at the pore entrance is 0.31 at t = 1 s and then increases to the peak value of 0.66 by t = 10 s; thereafter, particles are accumulated on the feed-side surface of the membrane and form the second and third cake layers on the membrane. The probability of particle existence inside the pore changes very little during the entire simulation,
and the average value between x = 6 and 10 at t = 15 s is 0.078. In summary, for ideal, numerically modeled straight cylindrical pores, the fouling of dm = 2.5d happens initially at the feed surface of the membrane without filling up in the pore. Thereafter, subsequent particles are accumulated sequentially, as shown in Fig. 8. The origin of this fouling is the feed surface of the membrane, and it then grows in the feed domain. As a result, the above explanations result that the number of particles sticking on the entrance and exit has great difference in the three kinds of pore size, as shown in Figs. 7(e), 8(e) and 13(b). However, there is no doubt that these situations will also depend on the inflow particle flux of the pore, that is, the particle concentration.
54
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
Fig. 9. Probability of particle existence for pore diameter dm = 3.6d.
Fig. 11. Permeate flux including a backwash operation for 1 s.
not a higher filling rate in the pores than that of dm = 3.6d, but the probabilities of particle existence are similar to that of dm = 3.6d (as shown in Figs. 9 and 14(a)). Therefore, the threshold pore size of the fouling mode exists between dm = 3d and 2.5d under this condition. 3.3. Simulation of the backwash operation Fig. 10. Probability of particle existence for pore diameter dm = 2.5d.
Considering the above, the correlation between pore size and fouling condition could be confirmed by the particle motion of two-way coupling model simulations under a fixed particle concentration (v = 5%). There are two fouling modes: in one, initially particles fill the pore and then form a cake layer on the surface of membrane; another is that particles are accumulated on the surface of the membrane without filling the pore and a cake layer grows across the entire filtration. The result of simulation of dm = 3d is
Simulation of a sequence of filtration including a backwash operation was performed after the normal filtration of dm = 3d for 30 s under a fixed particle concentration (v = 5%). The pressure of the backwash operation was eight times that of normal filtration, and the period of backwash operation was 1 s. Fig. 11 shows the permeate flux through a sequence of filtration including a backwash operation, and Fig. 12 shows snapshots of the backwash operation. In Fig. 12, the color of the particle shows its velocity. Just after beginning, t = 0.04 s, as is shown in Fig. 12(b), both particles in the pore and on the feed-side surface of the membrane have a finite velocity. The particles on the feed side of the membrane
Fig. 12. Snapshots of the backwash operation for 1 s.
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
form a large cluster and begin to move together simultaneously, as shown in Fig. 12(c); the particles in the pore follow. Then, as shown in Fig. 12(d), some particles remaining on the permeate-side surface of the membrane and in the pore form some small groups and grow away from the membrane. After the 1 s backwash operation,
55
some particles are left on both sides of the membrane, as shown in Fig. 12(f). Fig. 13 shows six snapshots of the start, middle, and end times before and after the backwash operation, and Fig. 14 shows the probabilities of particle existence. For the initial condition (t = 0 s)
Fig. 13. Snapshots of filtration before and after the backwash operation for dm = 3d: (a) 0 s, (b) 15 s, (c) 30 s, (d) 0 s after the backwash operation, (e) 15 s after the backwash operation, and (f) 30 s after the backwash operation.
56
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
Fig. 15. Normalized flux and total resistance as a function of time for filtration before and after the backwash operation.
Fig. 14. Probability of particle existence for filtrations before and after the backwash operation for dm = 3d: (a) first filtration before the backwash operation and (b) second filtration after the backwash operation.
of the first normal filtration, there are no particles attached to the membrane, as shown in Fig. 13(a). Conversely, for the initial condition of the second normal filtration just after the backwash operation, some particles are attached on the membrane, as shown in Fig. 13(d), and the probabilities of particle existence of both the entrance and the exit of the membrane are about the same value, 0.3, as shown in Fig. 14(b). In the case of idea, numerically modeled straight cylindrical pores, there are no particles left on the inside wall of the pore after the backwash operation. At the end of the first filtration, we can see that particles almost fill the entire pore, as shown in Fig. 13(c), but at the end of the second filtration, this
condition is not evident, as shown in Fig. 13(f). From these results, we observe that the fouling mode is different between the first and the second filtration. We can also confirm this difference from the normalized flux (J/J0 ) curve in Fig. 15, where J0 is the maximum value for each filtration. In addition, regarding the time-history of the probability of particle existence, Fig. 14(a) of the first filtration is similar to Fig. 9 of the dm = 3.6d filtration, and Fig. 14(b) of the second filtration is similar to Fig. 10 of the dm = 2.5d filtration. Additionally, the total resistance versus time of the second filtration draws a very similar curve to that of the dm = 2.5d filtration, as shown in Fig. 15. In fact, the detailed observation of short-interval snapshots shows that the second filtration of 3d is similar to the filtration of 2.5d in the conditions of the particle fouling. In the case of the membrane filtration process with straight cylindrical pores, the remaining particles on the feed-side surface of the membrane become the origin of the fouling after the backwash operation and form the cake layer on the feed-side surface of the membrane. This is supposed to be due to a decrease of the effective pore area. Indeed, it can be seen that the real area of the pore is reduced by the deposited particles, as shown in Fig. 16. We attempted to convert the binarized images of 0 and 1 from the grayscale images of Figs. 16(a) and (b) and analyze the total number of pixels in each value. The obtained ratio of the pore area (in black) to the total membrane area is 0.318. The effective pore area is assumed
Fig. 16. Attached particles on both sides of the membrane after the backwash operation: (a) feed side and (b) permeate side.
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
to be smaller than this value and closer to the porosity of pore size dm = 2.5d (ε = 0.306). The two-dimensional information (effective pore area) obtained from two faces of membrane after backwash operation is very useful because this study, which uses the ideal, numerically modeled pore, does not have the remaining particles in the pores. From the discussions in this study about the microfiltration with uniform pores including the after backwash operation, we speculate that the initial porosity εi (εi = (effective pore area)/(membrane area)) is a key factor to decide the fouling mode or condition. It is because new generated particles at the feed plane of 6d distance away from the membrane do not know the information of the membrane, such as the zeta potential and the number of sticking particles, in the initial stages of filtration and just move along a streamline with a speed according to the permeate flux. Nowadays, we cannot find out the related theory to permit direct comparison of our hypothesis. However, we think that the pore scale approach that Bacchin et al. [29] are recently working on shows great promise as a comparison with our results, and we are looking forward to the further study and advanced theory. 4. Conclusion We carried out numerical simulations of pressure-driven deadend microfiltration with a backwash operation using the two-way coupling model in consideration of particle–fluid interactions. This study numerically modeled membranes with regularly spaced straight pores, which are assumed to be track-etched pores. On the basis of the results obtained by microfiltration under a fixed particle concentration, the present paper examined the effect of pore size on particle fouling by comparing snapshots of particle motion, the permeate flux, and the total resistance of membranes over a period of time. Based on the results of simulations of a sequence of filtration including a backwash operation, we discussed the fouling conditions for filtration before and after the backwash operation. The results of microfiltration and the resulting particle motion confirm two modes for particle fouling: one is that initially particles fill the pore and form a cake layer on the surface of membrane, and another is that particles are accumulated on the surface of the membrane without filling the pore and a cake layer grows across the entire filtration. Here, based on a particular pore size, the larger pore size membrane shows the former mode, and the smaller pore size membrane shows the latter mode. Specifically, for the larger pore size membrane (dm = 3.6d), the total resistance versus time curve is initially concave-up due to pores fouling internally, and then the total resistance versus time curve is concave-down due to formation of the cake layer on the surface of membrane. For the smaller pore size membrane (dm = 2.5d), the total resistance versus time is concave-down for the entire simulation, and the fouling appears to be due to a growing cake layer on the surface of the membrane. These general behaviors support the hypothesis induced from experiments [27,28]. In addition, we compared the fouling conditions for filtration before and after the backwash operation. The particles remaining on the membrane after the backwash operation decrease the effective pore area and cause a change in fouling conditions from pore blockage to cake formation. In other words, in this study the backwash operation changed the fouling mode from the former mode of the larger pore size membrane to the latter mode of the smaller pore size membrane.
57
Acknowledgements The authors would like to thank Dr. Osamu Koike of the University of Tokyo for helpful discussions and comments. This work was performed as part of the Meso-Simulation Consortium. References [1] M. Cheryan, Ultrafiltration and Microfiltration Handbook, 2nd Eds., CRC Press, Boca Raton, Florida, 1998. [2] K.J. Howe, M.M. Clark, Fouling of microfiltration and ultrafiltration membrane by natural waters, Environ. Sci. Technol. 36 (2002) 3571–3576. [3] K. Akamatsu, W. Lu, T. Sugawara, S. Nakao, Development of a novel fouling suppression system in membrane bioreactor using an intermittent electric field, Water Res. 44 (2010) 825–830. [4] C.D. Orsello, W. Li, C.C. Ho, A three mechanism model to describe fouling of microfiltration membrane, J. Membr. Sci. 280 (2006) 856–866. [5] K.J. Hwang, C.Y. Liao, K.L. Tung, Effect of membrane pore size on the particle fouling in membrane filtration, Desalination 234 (2008) 16–23. [6] S. Kimura, S. Sourirajan, Analysis of date in reverse osmosis with porous cellulose acetate membrane use, AIChE J. 13 (1967) 497–503. [7] K.L. Tung, C.J. Chuang, Effect of pore morphology on fluid flow and particle deposition on a track-etched polycarbonate membrane, Desalination 146 (2002) 129–134. [8] M.M. Kim, A.L. Zydney, Effect of electrostatic, hydrodynamic, and Brownian forces on particle trajectories and sieving in normal flow filtration, J. Colloids Interface Sci. 269 (2004) 425–431. [9] M. Rahimi, S.S. Madaeni, M. Abolhasani, A.A. Alsairafi, CFD and experimental studies of fouling of a microfiltration membrane, Chem. Eng. Process: Process Intens. 48 (2009) 1405–1413. [10] M. Fujita, Y. Yamaguchi, Multiscale simulation method for self-organization of nanoparticles in dense suspension, J. Comput. Phys. 223 (2007) 108–120. [11] M. Fujita, Y. Yamaguchi, Simulation model of concentrated colloidal nanoparticulate flows, Phys. Rev. E 77 (2008) 026706. [12] T. Ando, K. Akamatsu, M. Fujita, S. Nakao, Direct simulation model of concentrated particulate flow in pressure-driven dead-end microfiltration, J. Chem. Eng. Jpn. 43 (2010) 815–828. [13] J. Hermia, Constant pressure blocking filtration laws-application to power law non-Newtonian fluid, Trans. Inst. Chem. Eng. 60 (1982) 183–187. [14] W.R. Bowen, J.I. Calvo, A. Hernandez, Steps of membrane blocking in flux decline during protein microfiltration, J. Membr. Sci. 101 (1995) 153–165. [15] L.D. Landau, E.M. Lifshitz, Fluid Mechanics, Pergamon Press, London, 1959. [16] K. Fujii, Unified zonal method based on the fortified solution algorithm, J. Comput. Phys. 118 (1995) 92–108. [17] T. Kajishima, S. Takiguchi, H. Hamasaki, Y. Miyake, Turbulence structure of particle-laden flow in vertical plane channel due to vortex shedding, JSME Int. J. Ser. B 44 (2001) 526–535. [18] Y. Nakayama, R. Yamamoto, Simulation method to resolve hydrodynamic interactions in colloidal dispersion, Phys. Rev. E 71 (2005) 036707. [19] O. Koike, S. Ohta, M. Fujita, Y. Yamaguch, Simulation model of concentrated colloidal rod-like nanoparticles, Jpn. J. Appl. Phys. 47 (2008) 8124–8130. [20] T. Tsuji, T. Tanaka, Ishida, Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe, Powder Technol. 71 (1992) 239–250. [21] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, McGraw-Hill, New York, 1970. [22] R.D. Minddlin, H. Deresiewicz, Elastic spheres in contact under varying oblique forces, J. Appl. Mech. 20 (1953) 327–344. [23] J.N. Israelachivili, Intermolecular and Surface Forces, Academic Press, London, 1992. [24] R.G. Cox, H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface – II. Small gap widths, including inertial effects, Chem. Eng. Sci. 22 (1967) 1753–1777. [25] A.J.C. Ladd, R. Verberg, Lattice-Boltzmann simulations of particle–fluid suspensions, J. Stat. Phys. 104 (2001) 1191–1251. [26] D.B. Spalding, A general purpose computer program for multi-dimensional oneand-two-phase flow, Math. Comput. Simul. 23 (1981) 267–276. [27] E.M. Tracy, R.H. Davis, Protein fouling of track-etched polycarbonate microfiltration membranes, J. Colloids Interface Sci. 167 (1994) 104–116. [28] C.C. Ho, A.L. Zydney, A combined pore blockage and cake filtration model for protein fouling during microfiltration, J. Colloids Interface Sci. 232 (2000) 389–399. [29] P. Bacchin, A. Marty, P. Duru, M. Meireles, P. Aimar, Colloidal surface interactions and membrane fouling: investigations at pore scale, Adv. Colloids Interface Sci. 164 (2011) 2–11.
Contents lists available at SciVerse ScienceDirect
Journal of Membrane Science journal homepage: www.elsevier.com/locate/memsci
Simulation of fouling and backwash dynamics in dead-end microfiltration: Effect of pore size Tsutomu Ando a,∗ , Kazuki Akamatsu b , Shin-ichi Nakao b , Masahiro Fujita c a b c
Department of Mechanical Engineering, College of Industrial Technology, Nihon University, 1-2-1 Izumi-cho, Narashino, Chiba 275-8575, Japan Department of Environmental and Energy Chemistry, Kogakuin University, 2665-1 Nakano-machi, Hachioji, Tokyo 192-0015, Japan Department of Chemical System Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
a r t i c l e
i n f o
Article history: Received 2 September 2011 Received in revised form 29 November 2011 Accepted 30 November 2011 Available online 9 December 2011 Keywords: Membrane pore size Particle fouling Dead-end microfiltration Backwash Simulation
a b s t r a c t We carried out numerical simulations of pressure-driven dead-end microfiltration with a backwash operation using a two-way coupling model with consideration of the particle–fluid interactions. This study numerically modeled a membrane with regularly spaced straight pores, which are assumed to be tracketched pores. On the basis of the results obtained by microfiltration under fixed particle concentration (5%), the present paper examines the effect of pore size on particle fouling by comparing snapshots of particle motion, the permeate flux, and the total resistance of the membrane over a period of time. The results of simulations confirm two modes for particle fouling when particle diameter d = 100 nm. The larger pore size membrane (dm = 3.6d) shows a fouling mode in which initially particles filled the pores and then form a cake layer on the surface of membrane. The smaller pore size membrane (dm = 2.5d) shows a fouling mode in which particles are accumulated on the surface of the membrane without filling the pores and a cake layer forms across the entire filtration. These general behaviors are in agreement with the hypothesis deduced from experiments. In addition, we carried out simulations of a sequence of filtration including a backwash operation and compared the fouling conditions for filtration before and after its inclusion. The particles remaining on the membrane after the backwash operation decrease the effective pore area and cause a change in fouling conditions from pore blockage to cake formation. © 2011 Elsevier B.V. All rights reserved.
1. Introduction In recent years, attention has been given to water treatment using membrane filtration, e.g., water purification, wastewater treatment and seawater desalination [1–3]. In general water treatment, nanofiltration and reverse osmosis primarily target the separation of ions, ultrafiltration targets the separation of molecules, and microfiltration targets the separation of particles and microorganisms. One of the most significant contemporary issues in microfiltration in the water treatment process is a decrease in membrane performance through fouling caused by particles suspended in water. As such, estimation of the decrement in microfiltration is essential in optimizing the membrane filtration process [4,5]. However, before arguing this issue, we should first and foremost understand the mechanism of fouling in microfiltration. What is the origin of fouling, and what are the precursors to said origin? There are many things we do not yet understand regarding fouling in microfiltration.
∗ Corresponding author. Tel.: +81 47 474 2338; fax: +81 47 474 2349. E-mail address: [email protected] (T. Ando). 0376-7388/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.memsci.2011.11.051
Since ultrafiltration and reverse osmosis are controlled by the diffusive motions of molecules, their performance can be partially estimated using the molecular diffusive model [6]. On the other hand, the hydrodynamic effect is dominant in microfiltration because particle motion is affected by solvent flow. Therefore, simulation of microfiltration requires the fluid dynamic model rather than the molecular diffusive model. Recently, many computational fluid dynamics (CFD) simulations have been used to evaluate the microfiltration process [7–9]. However, most of these extant CFD filtration simulations were based on the one-way method. This method first resolves the flow field around the membrane, and then particles are traced along the streamline in the flow field. In this case, particle motion does not affect solvent flow around the particles themselves; therefore, the one-way method cannot precisely resolve questions of concentrated particulate flow. The authors have worked toward developing a meso-scale model in which we consider the dynamics of Brownian particles and the electrostatic force between two particles. We have also conducted simulation studies of meso-scale concentrated particulate flow. This simulator, called “SNAP-F (Structure of NAnoParticles Flow),” is a two-way coupling model that can simultaneously resolve both the motions of particles and fluid while also incorporating particle–fluid interactions [10–12]. As the lattice spacing
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
49
used to resolve the flow field is much smaller than the diameter of particles and the flow of the surface on particles is resolved with high resolution, we can calculate the solvent flow through the fouling particles. SNAP-F applied to the microfiltration process can thereby estimate the filtration performance of a membrane system. Moreover, in addition to the process of surface deposition and the fouling of pores, SNAP-F can also simulate the backwash dynamics for the detachment of particles from both the surface of the membrane and the inside of the pores; this reveals the quantitative effect of the actual filtration process including both fouling and backwashing. Previously, in two-dimensional simulations of concentrated colloidal nanoparticles, we have shown that our model of solid contact interactions, including frictional forces between nanoparticles, describes the dynamics of aggregated colloidal nanoparticles [10], and we have discussed the aggregation and dispersion structure of nanoparticles (depending on the Peclet number) in simple shear flows [11]. In a previous paper [12], we developed a microfiltration membrane model with a single straight pore and carried out three-dimensional simulations of concentrated particulate flow in pressure-driven dead-end microfiltration. We discussed the filtration performance to show the permeate flux and particle rejection with visualization of particle fouling on a membrane model with a single straight pore. The present paper carries out simulations of pressure-driven dead-end microfiltration with a backwash operation. Here, we apply the membrane model with four regularly spaced cylindrical pores, which are assumed to be track-etched pores. This study addresses three kinds of membranes with more than twice and less than four-times pore size of the particle diameter. Since these pore diameters are larger than the particle diameter and are not large enough to accumulate on the inside wall of the pore, fouling is not classified by the conventional fouling model (such as “complete blocking” and “standard blocking”) [13,14]. On the basis of the results of microfiltration, we examine the effect of pore size on particle fouling by comparing snapshots of particle fouling, the permeate flux and the total resistance of membrane over a period of time. We then perform backwash simulation for the membrane with pore size three-times that of the particle diameter. Based on the results of the backwash operation, we discuss the fouling conditions for filtration before and after the backwash operation.
In the present paper, we used a two-step method [17,18]. The previous paper used a one-step method and fixed the pressure in an inner region of particles to zero in order to satisfy the conservation of mass on fluid–particle coupling, but the equation of continuity was not satisfied on the boundary of the particles. This study divides the pressure into two to satisfy the conservation of mass: the intermediate pressure without particles, p∗ , and the correcting pressure with particles, p .
2. Simulation model
∇ · vn+1 = 0
2.1. Flow of solvent
Then, operating on both sides of Eq. (6), and substituting it for Eq. (8), the Poisson equation of intermediate pressure is given by
The momentum equation of the entire flow field, containing both particles and solvent is the fluctuating Navier–Stokes equation [15], shown as follows [12].
∇ 2 p∗ =
∂v 1 1 1 + (v · ∇ )v = − ∇ p + ∇ 2 v + ∇ · S − D + ␣ ∂t ␣=
vp − v 1 1 + (v · ∇ )v − ∇ 2 v − ∇ · S + D t
p = p∗ + p
(3)
Substituting Eq. (3) for Eq. (1), we obtain the following equation using the flow velocity without particles v∗ .
1 1 1 ∂v∗ = (1 − ) −(v∗ · ∇ )v∗ + ∇ 2 v∗ + ∇ · S − D − ∇ p∗ ∂t +
vp − v∗ t
−
1 ∇ p
(4)
Here the acceleration of solvent received from the rigid body is defined by ␣ =
vp − v∗ t
+ (v∗ · ∇ )v∗ −
2 ∗ 1 1 1 ∇ v − ∇ · S + D − ∇ p (5)
Eq. (4) is divided into 2 steps predictor–corrector method as follows. (1st step)
v ∗ − vn t
= (1 − ) −(v∗ · ∇ )v∗ +
(2nd step)
vn+1 − v∗ = t
according
to
the
1 2 ∗ 1 1 ∇ v + ∇ · S − D − ∇ p∗
vp − v∗ t
−
1 ∇ p
(6)
(7)
For the 1st step, the continuity equation without particles is expressed by
∇ · v∗ = 0
(8)
For the 2nd step, the continuity equation of the new entire flow field with particles is expressed by
(9)
∇ · vn t
+∇ ·
(1 − )
−(v∗ · ∇ )v∗ +
2 ∗ 1 1 ∇ v + ∇ ·S− D
(10)
(1)
Then also, operating on both sides of Eq. (7), and substituting it for Eq. (9), the Poisson equation of correcting pressure is given by
(2)
∇ 2 p =
where t denotes time, v is the fluid velocity vector on the computational lattice, is the solvent density, p is the pressure, is the solvent viscosity, S is the thermal fluctuating stress tensor, and D is the pressure gradient vector. is the volume fraction of particles and membrane, and this value is [0,1] considering the overlap between two particles or between a particle and the membrane. ˛ is the acceleration vector associated with the rigid velocity of particles on the lattice. Our previous paper [12] defined the membrane using the FSA (fortified solution algorithm) [16], but the present paper treats the membrane the same as particles, i.e., the FSA is not used. Here the velocity of the membrane is fortified as zero.
∇ · v∗ t
vp − v∗
+∇ ·
t
(11)
The algorithm to evolve the entire flow field from n-step to n + 1step is given as follows: 1. Solve the prediction values v∗ , p∗ from the present value vn using Eqs. (6) and (10). 2. Solve the new value vn+1 and the correction value p from the prediction value v∗ using Eqs. (7) and (11). 3. Calculate the new value pn+1 from the prediction value p∗ and the correction value p using Eq. (3). 4. Calculate the force and torque acting on particles using ␣ of Eq. (5).
50
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57 Table 1 Simulation condition.
2.2. Motion of particles The motion of particles in solvent is governed by Newton’s equation of motion. The translational motion of the l-th particle is expressed by the following equation [12]. ml
dV l = F cl + F el + F vl + F hl + F lb l dt
Particle
(12)
where ml is the mass of the particle, Vl is the translational velocity vector of the particle, Fcl is the contact force, Fel is the electrostatic v
force, Fl is the van der Waals force,
Fhl
is the hydrodynamic force,
and Flb l is the lubrication force. The rotational motion of the l-th particle is expressed by the following equation. Il
dl = Tcl + Thl dt
(13)
where Il is the moment of inertia of the particle, l is the angular velocity vector of the particle, Tcl is the contact torque, and Thl is the hydrodynamic torque. The contact force and torque are calculated using the Voigt model and by considering the Hertzian theory and the Mindlin model [10,11,19–22]. Here, the normal and the tangential contact forces both between two particles and between a particle and the membrane are given by a combination of elastic force and damping force. Then, the electrostatic force and the van der Waals force both between two particles and between a particle and the membrane are given by the Derjaguin–Landau–Verwey–Overbeek (DLVO) theory [23]. In addition, as in the previous paper [12], the present study considers the normal component of the lubrication force between pairs of particles or between a particle and a plane wall when the gap width tends to zero [24,25]. Because the lubrication force is a dominant force on a scale that is smaller than the lattice resolution, the analytical expressions of the lubrication force both between particles and between a particle and the membrane have been added to compensate the underestimation by the lack of lattice resolution in this study. The hydrodynamic force F hl and torque Thl acting on the l-th particle are given as below, using the volume integration of ␣ .
F hl
l (r){␣ (r) + D}dV
(14)
l (r){r l (r) × ␣ (r)}dV
(15)
p
=− Vp
T hl = −
p
Vp
Solvent
where V p denotes the region occupied by the l-th particle. The diffusion of particle is obtained through the hydrodynamic force F h and torque T h including the stochastic fluctuating stress. This results in the many-body Brownian motion of the particles.
Membrane
Nomenclature
Value
Diameter,d [nm] Density,p [kg/m3 ] Young’s modulus,E p [Pa] Poisson ratio, [–] Coefficient of friction between particles, p [–] Coefficient of friction between particle and membrane, m [–] Zeta potential, p [mV] Concentration, v [%] Hamaker constant, A [J]
100 3.97 × 103 3.65 × 109 0.254 0.1
Density, [kg/m3 ] Viscosity, [Pa s] Relative permittivity, εr [–] Ionic strength of the electrolyte, I [mol/l] Valence of ions [–] Temperature, T [K]
1.0 × 103 1.0 × 10−3 80 10−4
Thickness, t m [nm] Pore diameter, dm [nm] Porosity, ε [–] Density, m [kg/m3 ] Young’s modulus, Em [Pa] Poisson ratio, m [–] Zeta potential, m [mV]
400 250, 300, 360 0.306, 0.441, 0.636 3.97 × 103 3.65 × 109 0.254 −40
0.1 −40 5 4.7 × 10−20
1 293.15
filtration region, as shown in Fig. 1(a), because we employ a periodic boundary condition on the faces parallel to the direction of flow. The computational domain is (x, y, z) = (16d, 8d, 8d), as shown in Fig. 1(b). The membrane is 4d (400 nm) in thickness with 4 regularly spaced straight cylindrical pores of standard diameter 3d (300 nm) and is located in the middle of the flow direction of the computational domain. This study has simulated the real microfiltration, but the membrane is thinner than the real size due to limitations of simulation time. The lattice spacing of the flow field is 1/9d. As the initial condition, there are some particles in the feed-side domain according to the set value of particle concentration. After the simulation begins, the flow is induced by a uniform pressure gradient, and particles are generated at a point in the y–z plane on the feed side using random numbers, depending on the particle concentration. At the boundaries on the feed and permeate side, which are perpendicular to the direction of flow, the Neumann condition ıvx /ıx = 0 is employed and vy = vz = 0. As for the time step, we set the value as t = 0.4 ns so as not to induce numerical instability in all cases. 3. Simulation results and discussion
2.3. Simulation method and conditions
3.1. Comparative simulation between 1 pore and 4 pores
The equation of translational and rotational motions of particles is solved by using the first-order Euler explicit scheme, and the equation of particle position is solved by using the Crank–Nicolson scheme. Then, the fluctuating Navier–Stokes equation is solved using the two-step method; the first step uses the Semi-implicit method for pressure-linked equations shortened (SIMPLEST) [26], and the second step uses the convergence algorithm by the iterative calculation of pressure. Table 1 shows the condition of filtration simulation in the present study, and Fig. 1 shows the computational domain. The solvent is assumed to be water, and the particle is a spherical body with a diameter d (d = 100 nm) with physical properties identical to those of an alumina particle. In the simulation of pressure-driven dead-end microfiltration, the computational domain is a part of the
Without particles (v = 0%), we carried out simulations of pressure-driven dead-end microfiltration of pore size dm = 3d (300 nm), which is an ideal track-etched pore. The pressure gradient Dx = −2 × 1011 Pa/m is applied due to limitation of simulation time, and the difference in pressure between the feed and the permeate surfaces of the membrane is about 320 kPa. Then, the solvent flow becomes the Hagen–Poiseuille flow in the pore and the maximum Reynolds numbers Rep (= d|vx |)/ at both S = 0 and S = / 0 are around 0.42 for both the 1-pore membrane and the 4-pore membrane under the same porosity (ε = 0.441). The flow accelerates just before the tracer particles reach the entrance of the pore and flow inward, decelerates rapidly at the pore exit, and diffuses radially in the permeate domain, as seen in the figure without particles in Ref. [12], where simulation of microfiltration
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
51
Fig. 2. Permeate flux as a function of time up to t = 10 s for the 1-pore membrane and the 4-pore membrane, pore diameter dm = 3d.
flux is given by the average flow rate in the pore divided by the cross-sectional area of membrane (8d*8d). As the permeate volume increases, the fluxes of all pore sizes become smaller and indistinguishable due to resistance from particle fouling becoming greater than the resistance of the membrane. The permeate flux data can be used to examine the fouling condition by replotting the data in terms of the total resistance Rt , as given by the following equation: Rt =
Fig. 1. Schematic of calculation domain for studying dead-end microfiltration: (a) periodic boundary condition on the parallel faces with respect to the direction of flow and (b) computational domain, including membrane, with four regularly spaced straight pores.
was performed using the membrane model with a single straight pore. Fig. 2 shows the permeate fluxes, J, of each pore obtained by two simulations of the 1-pore membrane and the 4-pore membrane of pore size dm = 3d (ε = 0.441), under a particle concentration of v = 5%. Here, the permeate flux J of each pore for the 4 pore is given by the average flow rate in the pore divided by the 1/4 cross-sectional area of membrane (4d*4d). As each pore for the 4pore membrane draws an independent flux curve, simulation with as many pores as possible is important to obtain the image of the whole membrane, as shown in Fig. 1(a). Therefore, this study uses the membrane model with 4 regularly spaced straight pores and discusses the effect of pore size on particle fouling in microfiltration.
P J
(16)
The results are shown in Fig. 4, where P is the transmembrane pressure (320 kPa). The total resistance for the smaller pore size membrane of dm = 2.5d is initially larger due to the smaller pores but only increases with a decreasing slope for the entire simulation. The total resistance of the larger pore size membrane of dm = 3d and 3.6d starts off lower but becomes larger than that of the smaller pore size membrane by the end of the microfiltration. The total resistance curves of the larger pore size membrane have a point of inflection and produce an S-shape. These general behaviors agree with the experimental results of Tracy and Davis [27].
3.2. Simulation of filtration for different pore sizes We carried out simulations of dead-end microfiltration for two kinds of pore size, dm = 3.6d (360 nm) and 2.5d (250 nm), under a fixed particle concentration of v = 5%. The porosities of each pore size membrane are ε = 0.636 and 0.306, and 1.5 and 1/1.5 times that of standard pore (dm = 3d, ε = 0.441), respectively. Regarding the three kinds of pore sizes, we compared the permeate flux versus total permeate volume, as shown in Fig. 3. Here the permeate
Fig. 3. Permeate flux as a function of permeate volume for three kinds of pore diameter dm .
52
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
Fig. 4. Total resistance as a function of time for three kinds of pore diameter dm .
We first observed snapshots of the fouling condition of the standard pore size dm = 3d, as shown in Fig. 5. The total resistance curve for dm = 3d increases with an increasing slope to the point of inflection (at around 5 s), after which resistance continues to increase but with a decreasing slope. Compared with Fig. 5(a)–(d), which are snapshots at 4, 5, 6, and 7 s, respectively, it can be seen that pores plug at 6 s. We can also confirm this fact from Fig. 6, which shows the corresponding probability of particle existence. t = 6 s is the time when the flux curve for dm = 3d intersects with that for dm = 2.5d in Fig. 3, that is, the flux J is 0.16 m3 /(m2 s). However, the fouling mode does not strictly alternate to either pore blockage or cake formation at t = 6 s, and, in fact, it is assumed that the cake formation starts from around t = 5 s of the inflection point. Next, we examined two simulations of dm = 3.6d as a lager pore and 2.5d as a smaller pore to compare the total resistance versus time curve. The snapshots of dm = 3.6d and 2.5d are shown in Figs. 7 and 8, respectively. For the larger pore size membrane (dm = 3.6d), the positive second derivative of the total resistance curve initially fouls pores internally, and then the negative second derivative of the total resistance curve forms the cake layer on the surface of the membrane, as shown in Figs. 4 and 7, respectively. This behavior is consistent with the discussions of Tracy and Davis [27] and Ho and Zydney [28]. For the smaller pore size membrane (dm = 2.5d), the second derivative of the total resistance curve is negative for the entire simulation, and the fouling appears to be due to a growing cake layer on the surface of the membrane, as shown in Figs. 4 and 8, respectively. This behavior is also consistent with the discussions of Tracy and Davis [27]. In addition, we quantitatively compare the time transition of the fouling condition to examine the corresponding probability of particle existence of dm = 3.6d and 2.5d, as shown in Figs. 9 and 10, respectively. As shown in Fig. 9, in the case of dm = 3.6d, by 1 s, the probability of particle existence between the entrance and the exit of the pore is already at nearly the same value, 0.22 and 0.23. On the other hand, at this time, the probability of particle existence inside the pore is less than half that of the entrance and the exit. At t = 2 s, the particle existence at the pore exit is a 79% increase from 1 s, but at the same time at the entrance, there is only a 15% increase, and that at the inside center of the pore (x = 8d) has changed very little. The probability of particle existence at t = 5 s draws a monotonically decreasing curve from the exit to the entrance of the pore. After t = 10 s, the probability of particle existence only increases in the feed side of the membrane. In summary, for ideal, numerically modeled straight cylindrical pores, the fouling of dm = 3.6d happens initially at the exit and thereafter
Fig. 5. Snapshots of simulation for pore diameter dm = 3d: (a) 4 s, (b) 5 s, (c) 6 s, and (d) 7 s.
propagates to the entrance and inside of the pore, as shown in Fig. 7. In other words, in the early stages of filtration without particles in the pore, the Poiseuille flow in the pore decelerates rapidly at the exit and diffuses radially in the permeate domain. Then, in tandem with this solvent flow, some rolling particles on the inside wall remain on the permeate-side surface of the membrane when
Fig. 6. Probability of particle existence for pore diameter dm = 3d.
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
53
Fig. 7. Snapshots of simulation for pore diameter dm = 3.6d: (a) 1 s, (b) 2 s, (c) 5 s, (d) 10 s, and (e) 15 s.
Fig. 8. Snapshots of simulation for pore diameter dm = 2.5d: (a) 1 s, (b) 2 s, (c) 5 s, (d) 10 s, and (e) 15 s.
particles vacate the pore. These remaining particles around the exit reduce the effective cross-sectional area of the pore and become the origin of fouling. As shown in Fig. 10, in the case of dm = 2.5d, at t = 1 s, the probability of particle existence at the pore exit is less than 0.1 and increases to 0.29 until t = 5 s; thereafter, it does not increase. Conversely, the probability of particle existence at the pore entrance is 0.31 at t = 1 s and then increases to the peak value of 0.66 by t = 10 s; thereafter, particles are accumulated on the feed-side surface of the membrane and form the second and third cake layers on the membrane. The probability of particle existence inside the pore changes very little during the entire simulation,
and the average value between x = 6 and 10 at t = 15 s is 0.078. In summary, for ideal, numerically modeled straight cylindrical pores, the fouling of dm = 2.5d happens initially at the feed surface of the membrane without filling up in the pore. Thereafter, subsequent particles are accumulated sequentially, as shown in Fig. 8. The origin of this fouling is the feed surface of the membrane, and it then grows in the feed domain. As a result, the above explanations result that the number of particles sticking on the entrance and exit has great difference in the three kinds of pore size, as shown in Figs. 7(e), 8(e) and 13(b). However, there is no doubt that these situations will also depend on the inflow particle flux of the pore, that is, the particle concentration.
54
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
Fig. 9. Probability of particle existence for pore diameter dm = 3.6d.
Fig. 11. Permeate flux including a backwash operation for 1 s.
not a higher filling rate in the pores than that of dm = 3.6d, but the probabilities of particle existence are similar to that of dm = 3.6d (as shown in Figs. 9 and 14(a)). Therefore, the threshold pore size of the fouling mode exists between dm = 3d and 2.5d under this condition. 3.3. Simulation of the backwash operation Fig. 10. Probability of particle existence for pore diameter dm = 2.5d.
Considering the above, the correlation between pore size and fouling condition could be confirmed by the particle motion of two-way coupling model simulations under a fixed particle concentration (v = 5%). There are two fouling modes: in one, initially particles fill the pore and then form a cake layer on the surface of membrane; another is that particles are accumulated on the surface of the membrane without filling the pore and a cake layer grows across the entire filtration. The result of simulation of dm = 3d is
Simulation of a sequence of filtration including a backwash operation was performed after the normal filtration of dm = 3d for 30 s under a fixed particle concentration (v = 5%). The pressure of the backwash operation was eight times that of normal filtration, and the period of backwash operation was 1 s. Fig. 11 shows the permeate flux through a sequence of filtration including a backwash operation, and Fig. 12 shows snapshots of the backwash operation. In Fig. 12, the color of the particle shows its velocity. Just after beginning, t = 0.04 s, as is shown in Fig. 12(b), both particles in the pore and on the feed-side surface of the membrane have a finite velocity. The particles on the feed side of the membrane
Fig. 12. Snapshots of the backwash operation for 1 s.
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
form a large cluster and begin to move together simultaneously, as shown in Fig. 12(c); the particles in the pore follow. Then, as shown in Fig. 12(d), some particles remaining on the permeate-side surface of the membrane and in the pore form some small groups and grow away from the membrane. After the 1 s backwash operation,
55
some particles are left on both sides of the membrane, as shown in Fig. 12(f). Fig. 13 shows six snapshots of the start, middle, and end times before and after the backwash operation, and Fig. 14 shows the probabilities of particle existence. For the initial condition (t = 0 s)
Fig. 13. Snapshots of filtration before and after the backwash operation for dm = 3d: (a) 0 s, (b) 15 s, (c) 30 s, (d) 0 s after the backwash operation, (e) 15 s after the backwash operation, and (f) 30 s after the backwash operation.
56
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
Fig. 15. Normalized flux and total resistance as a function of time for filtration before and after the backwash operation.
Fig. 14. Probability of particle existence for filtrations before and after the backwash operation for dm = 3d: (a) first filtration before the backwash operation and (b) second filtration after the backwash operation.
of the first normal filtration, there are no particles attached to the membrane, as shown in Fig. 13(a). Conversely, for the initial condition of the second normal filtration just after the backwash operation, some particles are attached on the membrane, as shown in Fig. 13(d), and the probabilities of particle existence of both the entrance and the exit of the membrane are about the same value, 0.3, as shown in Fig. 14(b). In the case of idea, numerically modeled straight cylindrical pores, there are no particles left on the inside wall of the pore after the backwash operation. At the end of the first filtration, we can see that particles almost fill the entire pore, as shown in Fig. 13(c), but at the end of the second filtration, this
condition is not evident, as shown in Fig. 13(f). From these results, we observe that the fouling mode is different between the first and the second filtration. We can also confirm this difference from the normalized flux (J/J0 ) curve in Fig. 15, where J0 is the maximum value for each filtration. In addition, regarding the time-history of the probability of particle existence, Fig. 14(a) of the first filtration is similar to Fig. 9 of the dm = 3.6d filtration, and Fig. 14(b) of the second filtration is similar to Fig. 10 of the dm = 2.5d filtration. Additionally, the total resistance versus time of the second filtration draws a very similar curve to that of the dm = 2.5d filtration, as shown in Fig. 15. In fact, the detailed observation of short-interval snapshots shows that the second filtration of 3d is similar to the filtration of 2.5d in the conditions of the particle fouling. In the case of the membrane filtration process with straight cylindrical pores, the remaining particles on the feed-side surface of the membrane become the origin of the fouling after the backwash operation and form the cake layer on the feed-side surface of the membrane. This is supposed to be due to a decrease of the effective pore area. Indeed, it can be seen that the real area of the pore is reduced by the deposited particles, as shown in Fig. 16. We attempted to convert the binarized images of 0 and 1 from the grayscale images of Figs. 16(a) and (b) and analyze the total number of pixels in each value. The obtained ratio of the pore area (in black) to the total membrane area is 0.318. The effective pore area is assumed
Fig. 16. Attached particles on both sides of the membrane after the backwash operation: (a) feed side and (b) permeate side.
T. Ando et al. / Journal of Membrane Science 392–393 (2012) 48–57
to be smaller than this value and closer to the porosity of pore size dm = 2.5d (ε = 0.306). The two-dimensional information (effective pore area) obtained from two faces of membrane after backwash operation is very useful because this study, which uses the ideal, numerically modeled pore, does not have the remaining particles in the pores. From the discussions in this study about the microfiltration with uniform pores including the after backwash operation, we speculate that the initial porosity εi (εi = (effective pore area)/(membrane area)) is a key factor to decide the fouling mode or condition. It is because new generated particles at the feed plane of 6d distance away from the membrane do not know the information of the membrane, such as the zeta potential and the number of sticking particles, in the initial stages of filtration and just move along a streamline with a speed according to the permeate flux. Nowadays, we cannot find out the related theory to permit direct comparison of our hypothesis. However, we think that the pore scale approach that Bacchin et al. [29] are recently working on shows great promise as a comparison with our results, and we are looking forward to the further study and advanced theory. 4. Conclusion We carried out numerical simulations of pressure-driven deadend microfiltration with a backwash operation using the two-way coupling model in consideration of particle–fluid interactions. This study numerically modeled membranes with regularly spaced straight pores, which are assumed to be track-etched pores. On the basis of the results obtained by microfiltration under a fixed particle concentration, the present paper examined the effect of pore size on particle fouling by comparing snapshots of particle motion, the permeate flux, and the total resistance of membranes over a period of time. Based on the results of simulations of a sequence of filtration including a backwash operation, we discussed the fouling conditions for filtration before and after the backwash operation. The results of microfiltration and the resulting particle motion confirm two modes for particle fouling: one is that initially particles fill the pore and form a cake layer on the surface of membrane, and another is that particles are accumulated on the surface of the membrane without filling the pore and a cake layer grows across the entire filtration. Here, based on a particular pore size, the larger pore size membrane shows the former mode, and the smaller pore size membrane shows the latter mode. Specifically, for the larger pore size membrane (dm = 3.6d), the total resistance versus time curve is initially concave-up due to pores fouling internally, and then the total resistance versus time curve is concave-down due to formation of the cake layer on the surface of membrane. For the smaller pore size membrane (dm = 2.5d), the total resistance versus time is concave-down for the entire simulation, and the fouling appears to be due to a growing cake layer on the surface of the membrane. These general behaviors support the hypothesis induced from experiments [27,28]. In addition, we compared the fouling conditions for filtration before and after the backwash operation. The particles remaining on the membrane after the backwash operation decrease the effective pore area and cause a change in fouling conditions from pore blockage to cake formation. In other words, in this study the backwash operation changed the fouling mode from the former mode of the larger pore size membrane to the latter mode of the smaller pore size membrane.
57
Acknowledgements The authors would like to thank Dr. Osamu Koike of the University of Tokyo for helpful discussions and comments. This work was performed as part of the Meso-Simulation Consortium. References [1] M. Cheryan, Ultrafiltration and Microfiltration Handbook, 2nd Eds., CRC Press, Boca Raton, Florida, 1998. [2] K.J. Howe, M.M. Clark, Fouling of microfiltration and ultrafiltration membrane by natural waters, Environ. Sci. Technol. 36 (2002) 3571–3576. [3] K. Akamatsu, W. Lu, T. Sugawara, S. Nakao, Development of a novel fouling suppression system in membrane bioreactor using an intermittent electric field, Water Res. 44 (2010) 825–830. [4] C.D. Orsello, W. Li, C.C. Ho, A three mechanism model to describe fouling of microfiltration membrane, J. Membr. Sci. 280 (2006) 856–866. [5] K.J. Hwang, C.Y. Liao, K.L. Tung, Effect of membrane pore size on the particle fouling in membrane filtration, Desalination 234 (2008) 16–23. [6] S. Kimura, S. Sourirajan, Analysis of date in reverse osmosis with porous cellulose acetate membrane use, AIChE J. 13 (1967) 497–503. [7] K.L. Tung, C.J. Chuang, Effect of pore morphology on fluid flow and particle deposition on a track-etched polycarbonate membrane, Desalination 146 (2002) 129–134. [8] M.M. Kim, A.L. Zydney, Effect of electrostatic, hydrodynamic, and Brownian forces on particle trajectories and sieving in normal flow filtration, J. Colloids Interface Sci. 269 (2004) 425–431. [9] M. Rahimi, S.S. Madaeni, M. Abolhasani, A.A. Alsairafi, CFD and experimental studies of fouling of a microfiltration membrane, Chem. Eng. Process: Process Intens. 48 (2009) 1405–1413. [10] M. Fujita, Y. Yamaguchi, Multiscale simulation method for self-organization of nanoparticles in dense suspension, J. Comput. Phys. 223 (2007) 108–120. [11] M. Fujita, Y. Yamaguchi, Simulation model of concentrated colloidal nanoparticulate flows, Phys. Rev. E 77 (2008) 026706. [12] T. Ando, K. Akamatsu, M. Fujita, S. Nakao, Direct simulation model of concentrated particulate flow in pressure-driven dead-end microfiltration, J. Chem. Eng. Jpn. 43 (2010) 815–828. [13] J. Hermia, Constant pressure blocking filtration laws-application to power law non-Newtonian fluid, Trans. Inst. Chem. Eng. 60 (1982) 183–187. [14] W.R. Bowen, J.I. Calvo, A. Hernandez, Steps of membrane blocking in flux decline during protein microfiltration, J. Membr. Sci. 101 (1995) 153–165. [15] L.D. Landau, E.M. Lifshitz, Fluid Mechanics, Pergamon Press, London, 1959. [16] K. Fujii, Unified zonal method based on the fortified solution algorithm, J. Comput. Phys. 118 (1995) 92–108. [17] T. Kajishima, S. Takiguchi, H. Hamasaki, Y. Miyake, Turbulence structure of particle-laden flow in vertical plane channel due to vortex shedding, JSME Int. J. Ser. B 44 (2001) 526–535. [18] Y. Nakayama, R. Yamamoto, Simulation method to resolve hydrodynamic interactions in colloidal dispersion, Phys. Rev. E 71 (2005) 036707. [19] O. Koike, S. Ohta, M. Fujita, Y. Yamaguch, Simulation model of concentrated colloidal rod-like nanoparticles, Jpn. J. Appl. Phys. 47 (2008) 8124–8130. [20] T. Tsuji, T. Tanaka, Ishida, Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe, Powder Technol. 71 (1992) 239–250. [21] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, McGraw-Hill, New York, 1970. [22] R.D. Minddlin, H. Deresiewicz, Elastic spheres in contact under varying oblique forces, J. Appl. Mech. 20 (1953) 327–344. [23] J.N. Israelachivili, Intermolecular and Surface Forces, Academic Press, London, 1992. [24] R.G. Cox, H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface – II. Small gap widths, including inertial effects, Chem. Eng. Sci. 22 (1967) 1753–1777. [25] A.J.C. Ladd, R. Verberg, Lattice-Boltzmann simulations of particle–fluid suspensions, J. Stat. Phys. 104 (2001) 1191–1251. [26] D.B. Spalding, A general purpose computer program for multi-dimensional oneand-two-phase flow, Math. Comput. Simul. 23 (1981) 267–276. [27] E.M. Tracy, R.H. Davis, Protein fouling of track-etched polycarbonate microfiltration membranes, J. Colloids Interface Sci. 167 (1994) 104–116. [28] C.C. Ho, A.L. Zydney, A combined pore blockage and cake filtration model for protein fouling during microfiltration, J. Colloids Interface Sci. 232 (2000) 389–399. [29] P. Bacchin, A. Marty, P. Duru, M. Meireles, P. Aimar, Colloidal surface interactions and membrane fouling: investigations at pore scale, Adv. Colloids Interface Sci. 164 (2011) 2–11.