Permeation of Dispersed Particles through a Pore and

Feb 25, 2013 - Center for Membrane and Film Technology, Department of Chemical Science and Engineering, Kobe University, 1-1 Rokkodai, Nada, Kobe ...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/IECR

Permeation of Dispersed Particles through a Pore and Transmembrane Pressure Behavior in Dead-End Constant-Flux Microfiltration by Two-Dimensional Direct Numerical Simulation Toru Ishigami,† Hiromi Fuse,† Shinichi Asao,‡ Daisuke Saeki,† Yoshikage Ohmukai,† Eiji Kamio,† and Hideto Matsuyama*,† †

Center for Membrane and Film Technology, Department of Chemical Science and Engineering, Kobe University, 1-1 Rokkodai, Nada, Kobe 657-8501, Japan ‡ Department of Mechanical Engineering, College of Industrial Technology, 1-27-1 Amagasaki, Hyogo 661-0047, Japan ABSTRACT: This paper presents numerical procedures for a coupled scheme that combines immersed boundary and discrete element methods to simulate two-dimensional fluid and particle motion. We simulated the dead-end, constant-flux microfiltration of dispersed particles and investigated the effect of hydrodynamic interactions between the fluid and particles on the transmembrane pressure. The permeation of a single particle through a pore was first simulated to obtain the fundamental characteristics of transmembrane pressure behavior. The effects of the particle dispersion volume fraction and the pore size on the transmembrane pressure were then investigated. The results show that the transmembrane pressure increases and fluctuates even in the absence of membrane fouling when a particle penetrates the pore. Transmembrane pressure behavior is also significantly affected by the number of particles within the pore and the location of particles around the pore. where Δp is the transmembrane pressure, Rfouling is the fouling resistance, J is the flux rate through the membrane, μf is the viscosity of the solvent, and Rm is the intrinsic membrane resistance, measured with pure solvent. Iritani et al.7 analyzed the transmembrane pressure behavior in dead-end microfiltration theoretically by applying the blocking law for constant-pressure filtration first proposed by Hermans and Bredee.8,9 The hydrodynamic interactions between the fluid and particles in the pores were not taken into account in these studies. The transmembrane pressure behavior may be influenced by the particle motions within the pore. However, to our knowledge, there have been no investigations on the relationship between hydrodynamic resistance and transmembrane pressure in dead-end constant-flux microfiltration. This is because it is very difficult to experimentally study the flow characteristics and particle motions within the pore that influence hydrodynamic interactions. Bungay and Brenner developed a theoretical model that represents increase in pressure associated with the motion of a spherical particle in a cylindrical pore.10 This model can be applied to a dilute particle system, but it cannot be applied to a concentrated particle system with particle−particle contact. In addition, they did not take into account the particle−particle and particle−pore surface physicochemical interactions. Numerical computations are effective for studying the flow field and particle motion inside the pore. A few researchers adopted numerical simulation to study the particle−fluid flow in microfiltration.11,12 Of particular interest, Ando et al. applied

1. INTRODUCTION Microfiltration (MF) plays an important role in water treatment processes. MF membranes are widely applied in drinking water purification, membrane bioreactors (MBRs), and pretreatment for seawater desalination because of their high rejection of particles and bacteria. Constant-flux filtration, by which end users can obtain a constant desired supply of filtered water, is often used in water treatment processes by MF. Defrance and Jaffrin1 reported the advantages of constantflux filtration for the performance of an MBR. In the case of the constant-flux filtration, the transmembrane pressure represents the water permeability of the MF membrane.1,2 From an energy-saving viewpoint, it is important in water treatment processes to understand the filtration dynamics of dispersed particles, because these govern transmembrane pressure behavior and this corresponds directly to the pumping power required to supply the feedwater. During dead-end MF, particles will be rejected but could also permeate the pores in some cases where the particle size is smaller than the membrane pore size. It is well-known that the transmembrane pressure increases through formation of a cake layer or adhesion of particles to the pore surface, known as membrane fouling.3 However, it is not as well understood that the transmembrane pressure also increases by the hydrodynamic interactions between fluid and particles, without fouling, as the suspension permeates through the membrane. Conventionally, the relationship between the transmembrane pressure and filtration resistance has been expressed by Darcy’s equation, as follows:2,4−6 J=

Received: Revised: Accepted: Published:

Δp μf (R fouling + R m)

(1) © 2013 American Chemical Society

4650

September 11, 2012 January 18, 2013 February 25, 2013 February 25, 2013 dx.doi.org/10.1021/ie302448x | Ind. Eng. Chem. Res. 2013, 52, 4650−4659

Industrial & Engineering Chemistry Research

Article

a direct numerical simulation technique to particulate flow in dead-end microfiltration.13,14 This simulation model, named Structure of NAno Particles (SNAP) and developed by Fujita et al.,15,16 is based on the discrete element method (DEM)17 for particle motions and computational fluid dynamics (CFD) for solvent flow, and takes into account the particle−particle and particle−-membrane physicochemical interactions and the hydrodynamic interactions between particle motions and fluid flow. This model was used to explore the filtration dynamics within a pore, but the authors dealt only with constant-pressure filtration, not constant-flux filtration. The aim of this study was to investigate the effects of fluid and particle motions on transmembrane pressure behavior during dead-end, constant-flux filtration. First, a direct numerical simulation model based on a two-way coupled scheme combining DEM and CFD was developed and validated. The effects of particle volume fraction in the feed solution and pore size on the transmembrane pressure were then investigated under constant-flux microfiltration. On the basis of these investigations, we clarify the effect of hydrodynamic interactions between the fluid and particles on the transmembrane pressure. This represents a first step toward achieving a fundamental understanding of the complex fluid dynamics of microfiltration.

momentum exchange is solved by the immersed boundary method18 without specification of the boundary condition around the particles. The volume-averaged velocity u is defined as follows: u = (1 − α)u f + α u p

where α is the local volume fraction of the particle and up is the velocity inside the particle. Here, the velocity can be decomposed using the translational and rotational components of the particle motion as follows: u p = vp + ωp × r (7) where r is the relative position vector between the center of the particle and the grid. The volume-averaged velocity is established with the assumption of nonslip at the interface. ∂u = −∇P + H + f ∂t

∂u f 1 + u f ·∇u f = − ∇p + νf ∇2 u f + f ∂t ρf

H = −u ·∇u + νf ∇2 u

Ip

dvp

u N + 1 = u N + Δt( −∇P + H + f)

dt

û = u N + Δt( −∇P + H) (3)

dωp dt

(4)

= Th + Tc

(5)

(10)

where N and Δt denote the time step and time increment, respectively. An approximate fluid velocity û in a cell is given as follows:

(2)

= Fh + Fe + Fv + Fc

(9)

A time advancement difference equation of eq 8 is described by the following equation:

(11)

If the cell is located inside the particle (α = 1), the velocity should be equal to the local particle velocity and f = (up − û)/ Δt. Conversely, if the cell is located within the fluid region (α = 0), the velocity should be equal to the local fluid velocity and f = 0. Thus, the relationship between the interaction force and local volume fraction can be expressed using a linear interpolation of α as follows:

where uf denotes the fluid velocity, t is the time, and p is the pressure. The density, ρf, and kinematic viscosity of the fluid, νf, were assumed to be constant. f is the interaction force between particle and fluid. A Cartesian grid system was used for the arrangement of the Eulerian variables. The motion of each particle was tracked with a Lagrange scheme by solving the equations of linear momentum and angular momentum as follows: mp

(8)

where P = p/ρf and H contains the convective and viscous terms as follows:

2. NUMERICAL METHOD 2.1. Governing Equations. The governing equations of fluid flow are the equation of continuity and the Navier−Stokes equation, as follows: ∇·u f = 0

(6)

f = α(u p − û )/Δt

(12)

The hydrodynamic force and torque acting on the particle are obtained by integrating the interaction force over the particle volume: Fh = −ρf

∫V f dV

(13)

∫V f × R dV

(14)

p

where mp denotes the mass of the particle, vp is the translating velocity of the particle, Ip is the inertia tensor, ωp is the angular velocity, and Fh, Fe, Fv, Fc, Th, and Tc are the hydrodynamic, electrostatic, van der Waals, and contact forces and the hydrodynamic and contact torques, respectively, acting on a particle, as described in the following sections. The Langevin equation is widely used in simulating particle motion, but the Brownian force is much smaller than the other interactions in the case of microfiltration, because the particle size is of micrometer order. Thus, the effect of the Brownian force was ignored in the present simulation. 2.2. Interaction Force between Particle and Fluid. The interaction force between the particle and fluid is described as the momentum exchange at the fluid−particle interface. The

Th = −ρf

p

where Vp is the volume of the particle and R is the relative position vector between the center and surface of the particle. 2.3. Particle Interaction. Contact forces are described by a Voigt model that consists of a spring, dashpot, and slider, based on DEM, as illustrated in Figure 1. Motions of all particles were tracked by the multibody interaction, based on the circle approximation.17 The membrane and particles can be assumed to overlap each other, and the contact force is calculated from the overlap distance between particles or particles and the membrane. The particle−particle or particle−membrane contact forces can be decomposed as follows: 4651

dx.doi.org/10.1021/ie302448x | Ind. Eng. Chem. Res. 2013, 52, 4650−4659

Industrial & Engineering Chemistry Research

Article

where r is the radius of the particle, A is the Hamaker constant, n is the number density of electrolyte ions, k is the Boltzmann constant, T is the temperature, H is the intersurface distance between particles, z is the electrolyte ionic valency, e is the elementary electric charge, ζ is the zeta potential, ε0 is the permittivity of a vacuum, and εr is the relative permittivity. 2.4. Numerical Method and Computational Conditions. The finite difference equivalents of the governing equations for fluid flow were solved numerically. The central finite difference and a hybrid difference scheme developed by Spalding21 were applied for diffusive and convective terms, respectively. A fully implicit form was used for a time marching method. The pressure field was calculated by the semi-implicit method for pressure-linked equations (SIMPLE) algorithm.22 The difference equations were solved iteratively, using the alternating direction implicit (ADI) method.23 The procedure for the particle motion was mathematically formulated by a Lagrangian approach, and a first-order explicit scheme was applied for the time-derivative term. Figure 2 shows the

Figure 1. Schematic of the DEM model for the numerical simulation of particles.

Fc n, t =

∑ (en,t + d n,t)

(15)

Fc s, t =

∑ (es, t + ds, t)

(16)

where e and d are the elastic force of the spring and the damping force of the dashpot, respectively. The subscripts “n” and “s” are normal and tangential components, respectively. The numerical procedure for calculating the contact force is based on the Hertz contact theory as en, t = en, t −Δt + K nΔun ,

Kn =

2bE , 3(1 − ν 2)

d n, t = ηn

ηn = 2 mK n

es, t = es, t −Δt + K sΔus ,

K s = K nc ,

Δun Δt

d s, t = ηs

Δus Δt

ηs = ηn c

(17)

(18)

(19)

Figure 2. Computational domain for the particulate flow in dead-end microfiltration.

(20)

where E is Young’s modulus, ν is the Poisson ratio, b is the length of contact line, c is the ratio of tangential spring constant to normal spring constant, and Δu is the relative displacement increment. The spring coefficient, K, and dashpot coefficient, η, are derived from Young’s modulus and the Poisson ratio. Derjaguin−Landau−Verwey−Overbeek (DLVO) theory19,20 considering the van der Waals and electrostatic interactions is applied to express the physicochemical interaction between particles or between a particle and the membrane as follows: ⎧ Ar − particle−particle ⎪ ⎪ 12H2 Fv = ⎨ ⎪− Ar particle−membrane ⎪ ⎩ 6H2

(21)

⎧ particle−particle ⎪ 2πrκεrε0ζ 2 ⎪ exp( −κH ) ⎪ ⎪ 1 + exp( −κH ) Fe = ⎨ ⎪ 2 particle−membrane ⎪ 4πrκεrε0ζ exp( −κH ) ⎪ ⎪ 1 + exp( −κH ) ⎩

(22)

computational domain for the simulation of dead-end, constant-flux filtration. Cartesian coordinates x and y were defined from the origin set at the bottom-left corner of the twodimensional domain. The size of the computational domain was Lx = 21.0d and Ly = 6.0d, where d is the diameter of 1.0 μm particles. The numbers of the main mesh were 281 × 81, and the time step size was 1.88 × 10−8 s. A membrane with a straight pore was mounted in the center of the streamwise direction of the computational domain, with the membrane thickness set at 10 μm. The thickness differed from the total thickness of a real membrane but did correspond to the skin layer thickness of an MF membrane.24 The fluid and particles flowed into the left boundary and exited from the right boundary. For boundary conditions, the inlet velocity on the feed side was uniform and constant under the assumption of constant-flux filtration. The Neumann condition was applied to the outlet boundary. The pressure at the outlet of the computational domain was fixed at p = 0. Steady velocity and pressure fields in the absence of particles were used as the initial conditions. The centers of particles were initially placed at x = 1.0 μm and then flowed in with the same translational velocity as the inlet uniform velocity of the fluid to satisfy the interparticle distance for the specified volume fraction. For the computational conditions, the Reynolds number based on the inlet velocity and the diameter of the particle, the zeta potential, and the Hamaker constant were set to 1, −40 mV, and 1.4 × 10−20 J, respectively. The physical properties of the particles, fluid, and membrane described by Ando et al.13 were used. Table 1 shows the parameters used for the calculation of the

where κ=

2nz 2e 2 εrε0kT

(23) 4652

dx.doi.org/10.1021/ie302448x | Ind. Eng. Chem. Res. 2013, 52, 4650−4659

Industrial & Engineering Chemistry Research

Article

shows the computational domain for the validation. The size of the computational domain was Lx = Ly = 15d. The particle was mounted at the center of the computational domain, and the calculation was continued until a steady state was reached. Figure 5 compares the drag coefficient obtained using the present method with the experimental results of Tritton et al.,25

Table 1. Parameters Used To Calculate Contact and DLVO Forces Young’s modulus [Pa] Poisson ratio [−] friction coefficient [−] zeta potential [mV] Hamaker constant [J] relative permittivity [−] electrolyte ionic valency [−] temperature [K]

3.83 × 109 0.34 0.1 −40 1.4 × 10−20 80 1 293.15

contact force and the DLVO force. The DLVO force curves between particles, calculated from the present conditions, are shown in Figure 3. The calculation result shows that

Figure 5. Comparison of the result obtained using the present method with the experimental result obtained by Tritton et al.25

showing that they are in good agreement. This indicates that the present method gives accurate solutions in the simulation. 3.2. Permeation Behavior of a Single Particle. To investigate the effect of particle motion on the velocity and pressure field within the pore, the constant-flux filtration of a single particle was simulated, with the pore size set to 2 μm. The time variation for the velocity contours is shown in Figure 6. The line on the particle is shown to help visualize the particle rotation. As shown in Figure 6a, no turbulence is observed in the initial flow field. In addition, Poiseuille flow forms inside the pore and the velocity gradient in the streamwise direction is sharp at the inlet and outlet of the pore. As the particle passes through the pore, the velocity around the particle decreases (Figure 6b−d). When the particle flows out of the pore, the velocity field recovers to the initial state, as shown in Figure 6e. Figure 7 shows the time variation of the pressure contours. When the particle enters the pore, the contour lines fluctuate from the particle surface (Figure 7b−d), indicating that the drag occurs between the particle and fluid. In addition, the pressure fluctuations are much larger when the particle is located at the inlet or outlet. These results indicate that the flow field is severely affected by the particle as it passes through the pore. Figure 8 shows the time variation of the transmembrane pressure. In this study, we define the transmembrane pressure as the pressure difference between the lower corner of the inlet and the outlet of the pore. It can be seen that Δp remains constant until around 2.2 μs and then shows peaks around 3.1 and 5.6 μs. The times of the first and second peaks correspond to those when the particle is located at the inlet and outlet of the pore, respectively. These tendencies are consistent with the pressure contours, as shown in Figure 7b,d. This is because a severe drag between the fluid and particle occurs when the particle passes across a region with a sharp velocity gradient. When the particle is located within the pore, the transmembrane pressure is slightly higher than that without the particle, as seen in the small pressure fluctuation in Figure 7c. However, the transmembrane pressure differed little from that

Figure 3. Relationship between van der Waals, electrostatic, and total potentials and the intersurface distance between particles under current conditions.

electrostatic repulsion is dominant if the intersurface distance exceeds around 2.0 × 10−2 μm but otherwise the van der Waals attraction is dominant in the present condition.

3. RESULTS AND DISCUSSION 3.1. Validation of Present Method. It is necessary to confirm the hydrodynamic interaction force between the particle and fluid, which is the primary force for transport of particles in microfiltration. To verify the interaction force between the fluid and the two-dimensional particle, the relationship between the Reynolds number Re, defined by the mean fluid velocity in the streamwise direction and the particle diameter, and the drag coefficient Cd was investigated. Figure 4

Figure 4. Computational domain for the flow around a single particle. 4653

dx.doi.org/10.1021/ie302448x | Ind. Eng. Chem. Res. 2013, 52, 4650−4659

Industrial & Engineering Chemistry Research

Article

Figure 7. Pressure contours for a single particle permeation, D = 2.0 μm.

Figure 6. Velocity contours for a single particle permeation, D = 2.0 μm.

without a particle when the particle was located outside the pore, because the particle does not significantly disturb the pressure field, as shown in Figure 7a,e. This transmembrane pressure increase due to particle motions shown in Figure 8 is much smaller than that in the later stage of membrane fouling in practical systems. However, in the initial stage of membrane fouling, the transmembrane pressure increases are small. Thus, the fouling mode in the initial stage of membrane fouling could be misjudged. The transmembane pressure increase due to particle motions within the pore is necessary for correctly understanding the fouling mode. Figure 9 shows the time variation of the drag force magnitude. The hydrodynamic drag force is much higher than the van der Waals and electrostatic forces shown in Figure 3. This indicates that the hydrodynamic interaction between the fluid and particle is dominant for particle motion in the present condition. In addition, the magnitude of hydrodynamic force increases until around 3.1 μs because the fluid velocity on the centerline of the channel increases as the fluid streamlines contract toward the pore. The subsequent hydrodynamic force decreases markedly and then remains almost constant as the particle enters the pore. Again, the drag force sharply increases and decreases slightly because of the sudden expansion of the channel at the outlet of the pore. On the basis of these investigations, we can conclude that transmembrane pressure behavior is significantly affected by the particle motion around the pore. 3.3. Effect of Volume Fraction on Transmembrane Pressure. The effects of the mean volume fraction of dispersed particles in the feedwater on the transmembrane pressure were

Figure 8. Time variation of transmembrane pressure for a single particle permeation, D = 2.0 μm.

investigated using the numerical method. The volume fraction in the feedwater, ϕ, was changed in three steps: 5.0 × 10−2, 1.0 × 10−1, and 1.5 × 10−1. The pore size of the membrane was fixed at 5.0 μm. Figure 10 shows the time variations of velocity contours. In all cases, no turbulence was observed during the initial states, as shown in Figure 10A1, B1, and C1. The velocity fields fluctuate around particles as they move toward the pore, while the velocity fields inside the pore barely change (Figure 10A2, B2, and C2). In addition, particles hardly rotate upstream of the pore. When the particles enter the pore, the flow field further fluctuates and the fluid velocity around the particles decreases (Figure 10A3, B3, and C3). Furthermore, it can be seen that 4654

dx.doi.org/10.1021/ie302448x | Ind. Eng. Chem. Res. 2013, 52, 4650−4659

Industrial & Engineering Chemistry Research

Article

feed side increases with time and increasing volume fraction, indicating that the transmembrane pressure increases. This is because the drag between the fluid and particles increases with increasing volume fraction, as described above. However, the pressure disturbance around the particles outside the pore is much lower than that inside the pore. This shows that the transmembrane pressure is significantly affected by particle motions inside the pore. The time variations in transmembrane pressure are shown in Figure 12. In all cases, the transmembrane pressures are similar to that without particles until around 2 μs. This is because the drag between the fluid and particles does not yet occur, as shown in Figure 11A1, B1, and C1. Thereafter, the transmembrane pressures show peak values when the first particles are positioned near the inlet of the pore. The values of peaks increase with time and reach a steady state. In addition, the transmembrane pressure increases with increasing volume fraction. These tendencies are consistent with the results of the velocity and pressure contours. The transmembrane pressure behaviors show periodicity in all cases. The peaks correspond to the times when particles are positioned at the inlet and outlet of the pore. Particles constantly enter from the inlet of the computational domain in the present conditions, giving rise to the transmembrane pressure periodicity. The lowest steady state values of transmembrane pressures increase with increasing particle volume fraction. Such increases in the transmembrane pressure are quite commonly seen when membrane fouling occurs.2,3,7 However, under the present conditions, particles do not block the pore and do not deposit on the pore surface and therefore no membrane fouling occurs. Even under the condition where there is no membrane fouling, the transmembrane pressure increases because of the drag between the fluid and particles10 and shows peak values when the first particles are positioned near the inlet of the pore, as shown in Figures 11 and 12. We believe that this is the first

Figure 9. Time variation in magnitude of the drag force acting on a single particle.

particles rotate when inside the pore because a higher shear rate acts on particles within the pore. When particles exit the pore, the flow field fluctuates downstream of the pore (Figure 10A4, B4, and C4). Compared with Figure 10A4, B4, and C4, the fluid velocity along the centerline of the pore decreases with increasing volume fraction, indicating that hydrodynamic resistance inside the pore increases. Figure 11 shows the time variations of pressure contours. No turbulence in the initial pressure fields was observed in any case (Figure 11A1, B1, and C1). However, the pressure field is significantly disturbed near the particle surfaces around the inlet of the pore (Figure 11A2, B2, and C2). When the particles enter the pore, pressure turbulence occurs continuously from the particle surfaces (Figure 11A3, B3, and C3). In addition, with increasing volume fraction, the pressure fluctuations around particles become more significant. This indicates that drag increases around the particle surfaces when the interparticle distance and the distance between particles and the pore surface become relatively small. The pressure on the

Figure 10. Velocity contours for different volume fractions in the feed solution, D = 5.0 μm. (A) ϕ = 5.0 × 10−2; (B) ϕ = 1.0 × 10−1; (C) ϕ = 1.5 × 10−1. 4655

dx.doi.org/10.1021/ie302448x | Ind. Eng. Chem. Res. 2013, 52, 4650−4659

Industrial & Engineering Chemistry Research

Article

Figure 11. Pressure contours for different volume fractions in the feed solution, D = 5.0 μm. (A) ϕ = 5.0 × 10−2; (B) ϕ = 1.0 × 10−1; (C) ϕ = 1.5 × 10−1.

gradient there. The pressure fields within the pore are then disturbed around the particles when they enter (Figure 14A3, B3, and C3) and then pass through the pore (Figure 14A4, B4, and C4). In addition, the pressures upstream of the pore increase slightly in all cases. Figure 15 shows the time variations of normalized transmembrane pressures, Δp/Δp0. Δp is normalized by the transmembrane pressure at steady state without particles (initial transmembrane pressure), because the absolute value of the initial transmembrane pressure varies significantly with different pore sizes under constant-flux-microfiltration conditions. To investigate the effect of particle motions on the transmembrane pressure, the effect of the initial pressure difference must be eliminated by normalizing the transmembrane pressure. The peak values in the case of D = 2.0 μm are much higher than those in other cases because the velocity gradient at the inlet and outlet of the pore becomes sharp with decreasing pore size, as shown in Figure 13. This sharp velocity gradient contributes a high drag force between the fluid and particles. However, the lowest values of Δp/Δp0 are similar. In all cases, the number of particles that permeate per unit time are similar because the volume fractions of feedwater and the water flux rates are the same in the present conditions. Thus, the effects of particles on the transmembrane pressure are similar, despite variations in pore size.

Figure 12. Time variations in transmembrane pressure with various volume fractions in the feed solution, D = 5.0 μm.

work that shows the effect of inlet and exit brings about the fluctuation of the transmembrane pressure. 3.4. Effect of Pore Size on Transmembrane Pressure. The effect of pore size on the transmembrane pressure was investigated. The pore size was changed in three steps of 2.0, 3.0, and 5.0 μm, and the volume fraction was set to 5.0 × 10−2. Figure 13 shows the time variations in the velocity contours. In all cases, particles move with the fluid flow and the fluid velocity around the particles then decreases when particles enter the pore. On the other hand, the central velocities in all cases do not decrease significantly as particles penetrate and pass through the pore, because of the low volume fraction of the feedwater. This result indicates that the effects of particle penetration on hydrodynamic resistance do not differ significantly with pore size. Figure 14 shows the time variations of pressure contours. In all cases, the pressure disturbance increases as particles approach the pore, as shown in Figure 14A2, B2, and C2. The pressure significantly fluctuates when the particles are located at the inlet of the pore because of the sharp velocity

4. CONCLUSIONS A direct numerical simulation model describing the motion of a particle dispersion in the pore was developed. This model is based on a two-way coupled scheme involving DEM and CFD that takes into account the hydrodynamic force between the fluid and particles, and the van der Waals, electrostatic, and contact forces between the particles and between the particles and the membrane. The present method was applied to constant-flux, dead-end microfiltration, and the relationship of the hydrodynamic force between the fluid and particles and the 4656

dx.doi.org/10.1021/ie302448x | Ind. Eng. Chem. Res. 2013, 52, 4650−4659

Industrial & Engineering Chemistry Research

Article

Figure 13. Velocity contours for different pore sizes, ϕ = 5.0 × 10−2. (A) D = 2.0 μm; (B) D = 3.0 μm; (C) D = 5.0 μm.

Figure 14. Pressure contours with different pore sizes in the case of ϕ = 5.0 × 10−2. (A) D = 2.0 μm; (B) D = 3.0 μm; (C) D = 5.0 μm.

transmembrane pressure was investigated. The following conclusions are drawn from the results. 1.The transmembrane pressure increases without membrane fouling when particles are positioned near the pore because of drag force between the fluid and particles. The transmembrane pressure shows peak values at the inlet and outlet of the pore, where the velocity gradients are relatively sharp. 2. As the volume fraction of the feedwater increases, the transmembrane pressure increases because of an increase in the number of particles in the pore. 3. A decrease in pore size caused an increase in the peak values of the transmembrane pressure because a sharp velocity gradient was formed at the inlet and outlet of the pore. On the

other hand, the lowest values were fairly similar because the number of particles within the pore per unit time were almost the same. It is important to note that a two-dimensional geometry such as those for the particles and pores in our simulations may overestimate the hydrodynamic forces. Thus, future work should extend these studies to a three-dimensional geometry. In addition, in this study, the Reynolds number, the particle size, the Hamaker constant, and the zeta potential were kept constant. These parameters may affect transmembrane pressure behavior. The effects of these parameters will also be investigated in the near future. 4657

dx.doi.org/10.1021/ie302448x | Ind. Eng. Chem. Res. 2013, 52, 4650−4659

Industrial & Engineering Chemistry Research

Article

Re = Reynolds number t = time, s T = temperature, K Tc = contact torque, N·m Th = hydrodynamic torque, N·m u = volume-averaged velocity vector, m·s−1 uf = fluid velocity vector, m·s−1 up = velocity vector inside particle, m·s−1 û = approximate fluid velocity vector, m·s−1 vp = translating velocity vector of particle, m·s−1 Vp = volume of particle, m3 x = horizontal coordinate, m y = vertical coordinate, m z = electrolyte ionic valency α = volume fraction of particle in the computational cell Δp = transmembrane pressure, Pa Δp/Δp0 = normalized transmembrane pressure Δt = time increment, s Δu = relative displacement increment, m ε0 = permittivity of vacuum, C·V−1·m−1 εr = relative permittivity η = dashpot coefficient, kg·s−1 ν = Poisson ratio νf = kinematic viscosity of fluid, m·s−2 ρf = density of fluid, kg·m−3 ωp = angular velocity vector, rad·s−1 ζ = zeta potential, V

Figure 15. Time variations in normalized transmembrane pressure with different pore sizes, ϕ = 5.0 × 10−2.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: +81-78-803-6180. Fax: +81-78-803-6180. Notes

The authors declare no competing financial interest.



NOMENCLATURE A = Hamaker constant, J b = length of contact line, m c = ratio of tangential spring constant to normal constant Cd = drag coefficient d = diameter of particle, m d = damping force vector, N D = diameter of particle, m e = elementary electric charge, C e = elastic force vector, N E = Young’s modulus, Pa f = hydrodynamic interaction vector between particle and fluid, m·s−2 Fc = contact force vector, N Fe = electrostatic force vector, N Fh = hydrodynamic force vector, N Fv = van der Waals force vector, N H = intersurface distance between particles, m H = convection and diffusion terms in Navier−Stokes equation, m·s−2 Ip = inertia tensor, kg·m2 J = flux rate, m3·m−2·s−1 k = Boltzmann constant, m2·kg·s−2·K−1 K = spring coefficient, N·m−1 Lx = horizontal length of computational domain, m Ly = vertical length of computational domain, m mp = mass of particle, kg n = number density of electrolyte ions p = pressure, Pa P = pressure divided by fluid density, m2·s−2 r = radius of particle, m r = relative position vector between center of particle and grid, m R = relative position vector between center and surface of particle, m Rfouling = fouling resistance, m−1 Rm = intrinsic resistance, m−1

Subscripts

c = contact e = electrostatic interaction f = fluid h = hydrodynamic interaction n = normal component p = particle v = van der Waals interaction s = tangential component Superscript



N = time step

REFERENCES

(1) Defrance, L.; Jaffrin, M. Y. Comparison between filtrations at fixed transmembrane pressure and fixed permeate flux: application to a membrane bioreactor used for wastewater treatment. J. Membr. Sci. 1999, 152, 203. (2) Ho, C.-C.; Zydney, A. L. Transmembrane pressure profiles during constant flux microfiltration of bovine serum albumin. J. Membr. Sci. 2002, 209, 363. (3) Palacio, L.; Ho, C.-C.; Zydney, A. L. Application of a PoreBlockage−Cake-Filtration Model to Protein Fouling During Microfiltration. Biotechnol. Bioeng. 2002, 79, 260. (4) Belfort, G.; Davis, H. D.; Zydney, A. L. The behavior of suspensions and macromolecular solutions in crossflow microfiltration. J. Membr. Sci. 1994, 96, 1. (5) Iritani, E. Modeling and evaluation of pore clogging of membrane in membrane filtration. Kagaku Kogaku Ronbunshu 2009, 35, 1. (6) Vyas, H. K.; Bennett, R. J.; Marshall, A. D. Influence of operating conditions on membrane fouling in crossflow microfiltration of particulate suspensions. Int. Dairy J. 2000, 10, 477. (7) Iritani, E.; Katagiri, N.; Sugiyama, Y. Clogging properties of membrane pores in constant-rate and constant-pressure microfiltration of dilute colloids. Kagaku Kogaku Ronbunshu 2011, 37, 323. (8) Hermans, P. H.; Bredée, H. L. Zur kenntnis der filtrationsgesetze. Recl. Trav. Chim. 1935, 54, 680.

4658

dx.doi.org/10.1021/ie302448x | Ind. Eng. Chem. Res. 2013, 52, 4650−4659

Industrial & Engineering Chemistry Research

Article

(9) Hermans, P. H.; Bredée, H. L. Principles of the mathematical treatment of constant-pressure filtration. J. Soc. Chem. Ind., London 1936, 55T, 1. (10) Bungay, P. M.; Brenner, H. Pressure drop due to the motion of a sphere near the wall bounding a Poiseuille flow. J. Fluid Mech. 1973, 60, 81. (11) Tung, K. L.; Chuang, C. L. Effect of pore morphology on fluid flow and particle deposition on a track-etched polycarbonate membrane. Desalination 2002, 146, 129. (12) Lu, W. M.; Tung, K. L.; Hwang, K. L. Effect of woven structure on transient characteristics of cake filtration. Chem. Eng. Sci. 1997, 52, 1743. (13) Ando, T.; Akamatsu, K.; Fujita, M.; Nakao, S. Direct simulation model of concentrated particulate flow in pressure-driven dead-end microfiltration. J. Chem. Eng. Jpn. 2010, 43, 815. (14) Ando, T.; Akamatsu, K.; Nakao, S.; Fujita, M. Simulation of fouling and backwash dynamics in dead-end microfiltration: Effect of pore size. J. Membr. Sci. 2012, 392−393, 48. (15) Fujita, M.; Yamaguchi, Y. Multiscale simulation method for selforganization of nanoparticles in dense suspension. J. Comput. Phys. 2007, 223, 108. (16) Fujita, M.; Yamaguchi, Y. Simulation model of concentrated colloidal nanoparticulate flows. Phys. Rev. E 2008, 77, 026706. (17) Cundall, P. A.; Strack, O. D. L. A discrete numerical model for granular assemblies. Geotechnique 1979, 29, 47. (18) Kajishima, T.; Takiguchi, S.; Hamasaki, H.; Miyake, Y. Turbulence structure of particle-laden flow in a vertical plane channel due to vortex shedding. JSME Int. J., Ser. B 2001, 44, 526. (19) Derjaguin, B. V.; Landau, L. D. Theory of the stability of strongly charged lyophobic sols and of the adhesion of strongly charged particles in solutions of electrolytes. Acta Physicochim. URSS 1941, 14, 633. (20) Verwey, E. J. W.; Overbeek, J. T. G. Theory of the Stability of Lyophobic Colloids; Elsevier: Amsterdam, 1948. (21) Spalding, D. B. A novel finite-difference formulation for differential expressions involving both first and second derivatives. Int. J. Numer. Methods Eng. 1972, 4, 551. (22) Patanker, S. V. Numerical Heat Transfer and Fluid Flow; Hemisphere Publishing: Washington, DC, 1980. (23) Peaceman, D. W.; Rachford, H. H. The numerical solution of parabolic and elliptic differential equations. J. Soc. Ind. Appl. Math. 1955, 3, 28. (24) Matsuyama, H.; Berghmans, S.; Lloyd, D. R. Formation of anisotropic membranes via thermally induced phase separation. Polymer 1999, 40, 2289. (25) Tritton, D. J. Experiments on the flow past a circular cylinder at low Reynolds numbers. J. Fluid Mech. 1959, 6, 547.

4659

dx.doi.org/10.1021/ie302448x | Ind. Eng. Chem. Res. 2013, 52, 4650−4659