NUMERICAL SIMULATION OF MICRO-FILTRATION OF OIL-IN-WATER EMULSIONS By Tohid Darvishzadeh A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering—Doctor of Philosophy 2014 ABSTRACT NUMERICAL SIMULATION OF MICRO-FILTRATION OF OIL-IN-WATER EMULSIONS By Tohid Darvishzadeh This study addresses the issue of oil removal from water using hydrophilic porous membranes. The effective separation of oil-in-water dispersions involves high flux of water through the membrane and, at the same time, high rejection rate of the oil phase. The effects of transmembrane pressure and crossflow velocity on rejection of oil droplets and thin oil films by pores of different cross-section are investigated numerically by solving the Navier-Stokes equation. We found that in the absence of crossflow, the critical transmembrane pressure, which is required for the oil droplet entry into a circular pore of a given surface hydrophilicity, agrees well with analytical predictions based on the Young-Laplace equation. An analytical expression for the critical pressure in terms of geometric parameters of the pore cross-section is validated via numerical simulations for a continuous oil film on elliptical and rectangular pores. With increasing crossflow velocity, the shape of the oil droplet is strongly deformed near the pore entrance and the critical pressure of permeation increases. We determined numerically the phase diagram for the droplet rejection, permeation, and breakup depending on the transmembrane pressure and shear rate. The critical pressure of permeation is identified as the line separating permeation and rejection regions. Using a novel method for computing the critical pressure, we investigated the effect of various physical and geometrical parameters on the critical pressure of permeation and breakup of droplets under shear flow. It is demonstrated numerically that the critical pressure of permeation increases with shear rate, viscosity ratio, surface tension coefficient, contact angle, and droplet size. On the other hand, droplet breakup at the pore entrance is facilitated at lower values of the surface tension coefficient, higher oil-to-water viscosity ratio, and larger droplet size. Using simple force and torque balance arguments, an estimate for the increase in critical pressure due to crossflow and the breakup capillary number is obtained and validated for different viscosity ratios, surface tension coefficients, contact angles, and drop-to-pore size ratios. ACKNOWLEDGMENTS Financial support from the Michigan State University Foundation Strategic Partnership Grant 71-1624 and the National Science Foundation (CBET-1033662) is gratefully acknowledged. Computational work in support of this research was performed at Michigan State University’s High Performance Computing Facility. iv TABLE OF CONTENTS LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 2 Details of numerical simulations . . . . . . . . . . . . . . . . . . . . . 10 Chapter 3 3.1 3.2 3.3 3.4 Effect of Crossflow Velocity and Transmembrane Pressure Microfiltration of Oil-in-Water Emulsions . . . . . . . . . . . . The Young-Laplace analysis for circular pores . . . . . . . . . . . . . . . . An oil droplet on a circular pore in the absence of crossflow . . . . . . . . Critical pressure for pores with arbitrary cross-section . . . . . . . . . . . 3.3.1 Thin oil film on the membrane surface with a rectangular pore . . 3.3.2 Thin oil film on the membrane surface with an elliptical pore . . . Sheared droplet on the membrane surface with a circular pore . . . . . . . 1 on . . . . . . . . . . . . . . . . . . . . . 22 22 24 26 29 32 34 Effect of Physical and Geometrical Parameters on Crossflow Microfiltration of Liquid-Liquid Emulsions . . . . . . . . . . . . . . Analytical formulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The effect of confinement on the critical pressure of permeation and breakup The effect of viscosity ratio on the critical transmembrane pressure . . . . . . The effect of surface tension on the critical pressure of permeation . . . . . . The effect of contact angle on permeation of the droplet . . . . . . . . . . . . The effect of droplet size on the critical pressure of permeation . . . . . . . . . . . . . . . 42 42 45 47 53 58 60 Chapter 4 4.1 4.2 4.3 4.4 4.5 4.6 Chapter 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 BIBLIOGRAPHY . . . . . . . . . . . . . v . . . . . . . . . . . . . . . . . . . 69 LIST OF FIGURES Figure 2.1 Values of the volume fraction for each of the phases and the interface (α = 1 in the first phase, α = 0 in the second phase, and 0 < α < 1 at the interface). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Figure 2.2 Profiles of a pinned droplet simulated using three levels of mesh resolution with the base level containing 30 grids along the diameter of the pore. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Figure 2.3 Deformation (left) and orientation angle of a droplet in simple shear flow for various capillary numbers: ● present simulations; ◇ VOF simulations of Li et. al; ◻ boundary integral simulations of Rallison; × boundary integral simulations of Kwak et. al; △ boundary integral simulations of Kennedy et. al; ○ experimental results of Rumschiedt and Mason; the solid line is for small deformation theory of Cox [78] (parameters: ρ1 = ρ2 = 1, µ1 = µ2 = 1, rd = 0.25, and γ˙ = 1) . . . . . . 16 Figure 2.4 Snapshots of the cross-section of a droplet in simple shear flow for various capillary numbers using present simulations (left) and the results of VOF simulations of Li et. al (right) [78] (parameters: ρ1 = ρ2 = 1, µ1 = µ2 = 1, rd = 0.25, and γ˙ = 1) . . . . . . . . . . . . . . . . . 17 Figure 2.5 Schematic of the hybrid mesh cross-section used in chapter 4. The mesh consists of coarse tetrahedral and fine hexagonal cells. . . . . . 18 Figure 2.6 Schematic representation of the oil droplet residing at the pore entrance in a rectangular channel with the corresponding boundary conditions. The width and length of the computational domain are fixed to 24 µm and 36 µm and the height of the domain is determined in Section 4.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 vi Figure 2.7 Schematic of the droplet cross-section on the pore. Critical pressure of permeation (P1 − P3 ) is calculated in three steps: (1) pressure jump across static interface (P1 − P2 ) is calculated from Young-Laplace equation, (2) pressure jump across the dynamic interface (P2 − P3 ) is computed numerically, and (3) pressure jump from steps 1 and 2 are added. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Figure 3.1 The critical pressure computed using Eq. (3.2) for the values of the pore radius 0.15 µm (continuous line), 0.2 µm (dashed line), and 0.3 µm (dotted line). Error bars are extracted from the numerical simulations (see text for details) with the parameters µo /µw = 2.4, ρo /ρw = 0.781, σ = 19.1 mN /m, and θ = 135○ . The symbols indicate (×) rejection and (○) permeation of the oil droplet. . . . . . . . . . . . . . . . . . . . 25 Figure 3.2 Schematic representation of the oil-water interface inside the pore of arbitrary cross-section with perimeter C and area A. The interface forms a constant angle θ with the inner surface of the pore. . . . . . 27 Figure 3.3 The critical permeation pressure for the oil film into the rectangular pore with different aspect ratios. The curves are Eq. (3.9) and the symbols are the numerical results for µo /µw = 2.4, ρo /ρw = 0.781, σ = 19.1 mN /m, and θ = 120○ . The symbols indicate (×) rejection and (○) permeation of the oil film. . . . . . . . . . . . . . . . . . . . . . 30 Figure 3.4 Snapshots of the oil-water interface inside the square pore for contact angles θ = 90○ , 120○ , 135○ , and 150○ . The critical pressure Eq. (3.9) is computed for the surface tension σ = 19.1 mN /m and the pore width 1.5 µm. Cases (a), (b), and (c) correspond to the stationary interface, while in the case (d) the interface is in transient state (see text for details). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Figure 3.5 The critical pressure of permeation for the oil film into the elliptical pore as a function of the major and minor radii. The curves are computed using Eq. (3.11). The symbols represent numerical results for the parameters µo /µw = 2.4, ρo /ρw = 0.781, σ = 19.1 mN /m, and θ = 120○ . The symbols indicate (×) rejection and (○) permeation of the oil film. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 vii Figure 3.6 Schematic representation of the oil droplet in the channel with the circular pore. The shear flow is induced by the upper wall moving with a constant velocity parallel to the stationary lower wall. The droplet radius is rd = 0.9 µm and the pore radius is rp = 0.2 µm. Periodic boundary conditions are applied at the inlet and outlet, while a constant pressure is maintained at the side walls. . . . . . . . . . . . 35 Figure 3.7 The phase diagram of the transmembrane pressure versus shear rate for the oil droplet with rd = 0.9 µm and the circular pore with rp = 0.2 µm. The contact angle is θ = 135○ . Each symbol represents a separate simulation that corresponds to either (◯) permeation, (square) breakup or (∇) rejection. Letters A, B, C, D, and E indicate operating conditions for the series of snapshots shown in Fig. 3.8. . . . . . . 36 Figure 3.8 Sequences of snapshots of the sheared droplet on the pressurized pore for different operating conditions as indicated in Fig. 3.7. The droplet radius is rd = 0.9 µm, the pore radius is rp = 0.2 µm, and the contact angle is θ = 135○ . The letters denote (A) rejection at low shear rates (∆t ≈ 15 µs), (B) permeation (∆t ≈ 15 µs), (C) rejection at high shear rates (∆t ≈ 2 µs), (D) local breakup (∆t ≈ 2.5 µs), and (E) breakup with necking (∆t ≈ 4 µs). . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Figure 3.9 The leakage volume as a function of the applied pressure for rd = 0.9 µm, rp = 0.2 µm, and γ˙ = 5 × 105 s−1 . The square symbols indicate the numerical results and the dashed line is the best fit to the data. The inset shows a snapshot of the process shortly after breakup of the droplet. The contact angle of the oil droplet in water is θ = 135○ . 39 Figure 4.1 The cross-sectional profiles of oil droplets in steady shear flow for the indicated confinement ratios when the capillary number is Ca = 0.021. The droplet radius is rd = 2 µm, the pore radius is rp = 0.5 µm, the contact angle is θ = 135○ , the surface tension coefficient is σ = 19.1 mN /m, and the viscosity ratio is λ = 1. . . . . . . . . . . . . . . . 46 Figure 4.2 The critical (breakup) capillary number as a function of the confinement ratio Hd /Hch . Other parameters are the same as in Fig. 4.1. . . 47 viii Figure 4.3 The percent increase in critical pressure of permeation as a function of the capillary number Ca = µw γr ˙ d /σ for the indicated viscosity ratios λ = µo /µw . Typical error bars are shown on selected data points. For each value of λ, the data are reported up to the critical capillary number above which droplets break into two segments. The droplet and pore radii are rd = 2 µm and rp = 0.5 µm, respectively. The contact angle is θ = 135○ and the surface tension coefficient is σ = 19.1 mN /m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 Figure 4.4 The cross-sectional profiles of the oil droplet residing on the circular pore with rp = 0.5 µm for the indicated capillary numbers. The viscosity ratio is λ = 1. Other parameters are the same as in Fig. 4.3. 50 Figure 4.5 The breakup time versus deformation time scale µw rd (1 + λ)/σ for the tabulated values of the viscosity ratio λ = µo /µw . Other system parameters are the same as in Fig. 4.3. The straight line is the best fit to the data. The error bars for the breakup time are about the symbol size. The inset shows the droplet profiles just before breakup for the same viscosity ratios. . . . . . . . . . . . . . . . . . . . . . . . . 51 Figure 4.6 The normalized percent increase in critical pressure of permeation versus the modified capillary number Ca fD (λ) for the selected values of the viscosity ratio λ = µo /µw . The functions fD (λ) and fT (λ) are given by Eq. (4.3) and Eq. (4.6), respectively. . . . . . . . . . . . . . . 52 Figure 4.7 The critical pressure of permeation as a function of shear rate for the indicated surface tension coefficients. The symbols (×) denote the analytical predictions of Eq. (3.2). The droplet and pore radii are rd = 2 µm and rp = 0.5 µm, respectively. The viscosity ratio is λ = 1 and the contact angle is θ = 135○ . . . . . . . . . . . . . . . . . . . . . . 54 Figure 4.8 The cross-sectional profiles of the oil droplet above the circular pore for the listed values of the surface tension coefficient. In all cases, the shear rate is γ˙ = 1.5 × 105 s−1 . Other parameters are the same as in Fig. 4.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 ix Figure 4.9 The breakup time versus deformation time scale µw rd (1+λ)/σ for the surface tension coefficients in the range from 9.55 mN /m to 38.2 mN /m. Other parameters are the same as in Fig. 4.7. The straight solid line is the best fit to the data. The cross-sectional profiles of the oil droplet just before breakup are displayed in the inset. . . . . . . . . . . . . . . 56 Figure 4.10 The percent increase in critical pressure of permeation as a function of the capillary number Ca = µw γr ˙ d /σ for the selected values of the surface tension coefficient. The rest of the material parameters are the same as in Fig. 4.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Figure 4.11 The critical pressure of permeation as a function of the capillary number for the indicated contact angles. The critical pressure at zero shear rate, given by Eq. (3.2), is denoted by the symbols (×). The droplet radius, pore radius, surface tension coefficient, and viscosity ratio are rd = 2 µm, rp = 0.5 µm, σ = 19.1 mN /m, and λ = 1, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Figure 4.12 The cross-sectional profiles of the oil droplet above the circular pore for the listed values of the contact angle when Ca = 0.022. Other parameters are the same as in Fig. 4.11. . . . . . . . . . . . . . . . . . 60 Figure 4.13 The critical pressure of permeation as a function of shear rate for the selected drop-to-pore size ratios. The symbols (×) indicate the critical pressure in the absence of flow calculated from Eq. (3.2). The pore radius, surface tension coefficient, contact angle, and viscosity ratio are rp = 0.5 µm, σ = 19.1 mN /m, θ = 135○ and λ = 1, respectively. 61 Figure 4.14 The difference in the critical pressure, Pcr − Pcr0 , versus the modified capillary number for five drop-to-pore size ratios r = rd /rp and rp = 0.5 µm. Other parameters are the same as in Fig. 4.13. The crosssectional profiles of the droplet just before breakup are shown in the inset for same r. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 x Chapter 1 Introduction With the recent advances in environmental and biological technologies, there has been increasing interest in characterization and modeling flows at the micron scales including flows in microchannels and nanochannels [1, 2], multiphase flows through porous media [3, 4], and droplet-based microfluidics [5, 6]. The industrial applications include oil extraction from porous media [7, 8, 9], treatment of oily wastewater [10, 11, 12], and encapsulation of molecules, cells, and microorganisms [13, 14, 15]. In recent years, microfluidics has emerged as a promising new tool to address various problems in different areas of research and industry including biology, petroleum industry, and environmental engineering [1, 16, 17]. Microfluidic devices are used as efficient analysis systems for DNA characterization and cell sorting processes in biomedical applications [18]. Moreover, measurements of chemical composition and properties of oil produced during drilling and production has enhanced significantly by using microfluidic sensors [19]. On the other hand, the science of microfluidic behavior of fluids and drops has contributed to the design of micro-separation systems to treat oily wastewater produced in oil spills, oil leaks, and industrial discharge [20]. In most of these processes, one or more phases are dispersed in a continuous phase in the form of emulsions, which are usually produced by shearing two immiscible phases against each other in the presence of surfactants [21]. Emulsions are used in a handful of different applications ranging from petroleum industry [22] to biomedical treatments [23]. In some cases, emulsions serve as means of transport of molecules, bio-reagents, and drugs, and ultimately provide the envi1 ronment for enhanced reactions [14, 24]. Another technological application of emulsions is to improve the transportability or displacement of highly viscous liquids. For example, heavy crude oil is emulsified to form a less viscous mixture to facilitate its transportation [25, 26]. In addition, oil-in-water emulsions are used to enhance recovery and increase sweep efficiency from crude oil reservoirs by blocking highly permeable paths and preventing channeling of the displacing fluid [7]. A very interesting part of microfluidics deals with the formation [27], transport [14], and stability [28] of emulsions. Depending on applications, producing stable emulsions could be either beneficial or problematic [29]. In food industries, it is essential for the emulsion to remain stable to retain the required flavor and quality [30]. On the other hand, emulsions with high stability produced during bio-diesel washing, crude oil extraction, and industrial wastewater are difficult to demulsify [29]. Common methods for separation of emulsions include evaporation of the continuous phase [31], destruction (demulsification) [32], and membrane filtration [10]. Micro-porous membranes are often used as an efficient tool to both produce emulsions [33] and to separate them from their bulk fluid [10]. The typical pore size for both applications is in the range of 0.1 µm to 10 µm [34]. Apart from conventional methods such as ultrasound emulsification [35] and stirring vessels [36], emulsions are produced through a process known as membrane emulsification [37]. Membrane emulsification generally involves a process in which micon-sized droplets are produced by forcing a liquid stream through a membrane pore into a channel in which another liquid is flowing. The emerging droplet breaks as a result of interactions between the shear, pressure, and capillary forces [38]. Membrane emulsification is more cost efficient, requires less energy density, and ultimately produces a narrower droplet size distribution than conventional methods [39, 40]. On the other hand, membrane microfiltration has proven to be an efficient way for separating oil-in-water emul2 sions [41, 42]. In comparison with the conventional methods of filtration (gravity separators, centrifuges, etc.), membrane microfiltration has several distinct advantages including reduced space requirements, higher permeate quality, and lower operating costs [10]. Despite its advantages, microfiltration efficiency can be greatly reduced because of membrane fouling at highly concentrated emulsions or long filtration times [43]. Fouling is generally caused by the accumulation of the rejected phase on the surface of the membrane or inside the pore. There are four main mechanisms (blocking laws) for membrane fouling, i.e., standard blocking, complete blocking, intermediate blocking, and cake formation [44]. Complete blocking is common for very dilute mixtures and during the initial stages of filtration when some pores are sealed by droplets and particles, thus reducing the permeate flux [45, 46]. Accumulation of the rejected droplets on the membrane surface results in the formation of the so-called cake layer, which is sometimes referred to as the secondary membrane as it adds a hydraulic resistance to the microfiltration process [47, 48, 49]. This mechanism is dominant at the final stages of filtration when the water flux depends mainly on the thickness of the cake layer. The present study is mainly related to the “complete blocking” mechanism in which a droplet blocks a pore completely and prevents the flux of the continuous phase [50]. The efficiency of the microfiltration process is determined by the properties of the membrane material and oil-in-water mixtures. For example, the permeate flux is highly dependent on the oil concentration, stability of the oil phase in water, and the size distribution of oil droplets [10, 51]. Moreover, the membrane properties such as membrane material, pore size and morphology, and membrane geometry affect the permeate flow resistance [10, 20]. It was shown that slotted (rectangular) pores resulted in higher flux rates compared to circular pores for similar operating conditions because of the lower fowling rate of the slotted-pore membranes [52]. Another approach to reduce the fouling rate is to introduce crossflow above 3 the membrane surface. This method, known as the “crossflow microfiltration”, reduces fouling by sweeping away the deposited drops and particles and, hence, decreases the thickness of the cake layer. Therefore, crossflow microfiltration systems tend to produce higher permeate fluxes for longer times compared to dead-end microfiltration systems [53, 54]. One of the aims of the present study is to investigate numerically the entry dynamics of oil droplets into a membrane pore in the presence of crossflow. Depending on the application and droplet size distribution, crossflowing systems could be either confined or unconfined. Most crossflow microfiltration processes use unconfined systems to produce a higher flux of the continuous phase. On the other hand, emulsification processes often use confined systems (e.g. t-junctions) to have a better control on the size and pattern of the emulsions generated [55]. During the last decade, a number of studies have investigated the process of droplet formation using cross-flowing streams in T-shaped junctions [56, 57, 58, 59]. In these microfluidic systems, two immiscible liquids are driven through separate channels until their streams meet at a junction, where the dispersed liquid extends into the continuous stream, resulting in periodic formation of equal-sized droplets [40]. Regardless of the specific channel geometry and wettability of the channel walls, breakup of the emerging droplet in a cross-flowing stream is determined by the viscous drag when the droplet remains unconfined by the microchannel [40]. For unconfined T-junctions, it was demonstrated experimentally that the droplet size strongly depends of the crossflow rate of the continuous phase and only weakly on the flow rate of the dispersed phase [57, 58]. It was also shown that, for given value of continuous phase flow rate, the size of oil droplets decreases with increasing viscosity ratio of the oil and water [58]. For microsized systems, the main force acting on a droplet in a confined environment is the upstream pressure, whereas for unconfined systems, the main contribution is from the shear stress on 4 the droplet. Consequently, breakup in unconfined systems is mainly shear dominated. However, droplets in confined channels, break as a result of blockage pressure gradient produced in the channel [40, 60]. Studying emulsions in general has been performed in two different scales, i.e. macroscopic or bulk scales and mesoscopic scales [22]. Early works in membrane emulsification and microfiltration includes performing bulk experiments to calculate average quantities and formulate empirical relations [33]. Bulk emulsification tends to produce broad distribution of droplet sizes as the shear stresses are not homogeneous [55, 14]. In recent years, however, with the development of advanced imaging and measuring techniques and computational resources, researchers studied a single droplet in a flow and related its microscopic behavior to the bulk properties of the flow [61, 62]. Macroscopic studies consider parameters such as droplet size distribution, dispersed phase concentration, and bulk properties such as permeate flux [63, 64]. On the other hand, the microscopic behavior of a droplet is quantified using different parameters like deformation, orientation angle, and breakup criteria [62]. In microfiltration of emulsions, there has been a limited number of studies investigating the connection between the macroscopic and mesoscopic variables [65, 66]. Micro-scale study of single droplets under various flow conditions has gained considerable attention in the past century. Taylor was the first to systematically study the deformation and breakup of a single droplet under shear flow. He performed experiments and developed theories explaining effect of viscosity on the droplet behavior [67, 68]. After Taylor, several other researchers have contributed theoretically [69, 70] and experimentally [71, 72] to the characterization of droplet deformation and breakup under shear flow. Bentley et. al [73] used a computer-controlled four-roll mill to investigate the effect of shear and extensional flows on a single droplet. Many researchers also performed numerical simulations using 5 boundary integral method to analyze droplets under different types of flow [62, 74]. Boundary integral methods have the advantage that the flow is solved only on the boundary of the droplet and the surface rather than the whole flow field and so the interface can be very accurately modeled. [75]. With the development of powerful computers, new algorithms were developed for flow solution and interface tracking with very high resolutions. One of these effective algorithms is the use of Volume of Fluid method to define the interface. This method tracks the interface between the two phases using a color function [76, 77]. The Volume of Fluids method has been extensively used to validate results of experiments and boundary integral simulations and to predict flow conditions [78, 79]. The deformation and orientation angle of droplets deposited on a surface under shear flow is greatly influenced by its contact angle, viscosity ratio, and the contact angle hysteresis [80]. Studying effect of shear flow on particles and droplets near walls has been received considerable attention in the past 50 years. Saffman concluded that a droplet near wall under shear flow could experience an upward lift force if inertial effects are not negligible [81]. O’Neill derived an exact solution for the linear Stokes flow over a spherical particle on a solid surface [82]. Price found the drag force acting on a hemispherical bump as a result of linear shear stokes flow on a solid wall [83]. Pozrikidis extended Price’s work to a spherical bump with arbitrary angle using boundary integral method [84]. Sugiyama et al. extended the work of Price by varying the viscosity ratio to values other that infinity (as considered by Price [83]) and found an exact solution for the linear shear flow past a hemispherical droplet on a solid surface. Assuming that the droplet is pinned to the surface, they computed the drag force, torque, and the deformation angle as function of the viscosity ratio [85]. The dynamics of droplet breakup in steady shear flow is determined by the relative competition of the viscous stress, pressure, and interfacial tension [62]. In general, the 6 breakup process is initiated by the droplet deformation, which is linearly proportional to the rate of shear [72]. When the critical deformation is reached, the droplet assumes an unstable configuration and undergoes a transient elongation before it breaks up [62]. It was also shown that the geometric confinement as well as the viscosity ratio of the dispersed and continuous phases influence droplet breakup [86]. In recent years, the problem of droplet deformation and breakup has been extensively studied numerically using Lattice Boltzmann [56, 87, 88], boundary integral [89, 90], and Volume of Fluid (VOF) [77, 78, 91] methods. The VOF method, used in the present study, has proven to be a powerful and efficient interface tracking algorithm that is both conceptually simple and relatively accurate [92]. Due to the conservative discretization of the governing equations in the VOF method, the mass of each fluid is accurately conserved [78, 93]. Furthermore, the ability of the VOF method to automatically capture local and global changes of the interface topology, e.g., coalescence and breakup of droplets, has made it attractive for various multiphase flow applications [92]. When studying micron-sized droplet deformation and breakup, there are multiple parameters to consider that are related to each other through dimensionless numbers. These parameters include viscosity of the continuous phase, surface tension coefficient between the two fluids, contact angle between the fluid-fluid interface and a solid surface attached to them, the viscosity ratio of the dispersed phase and the continuous phase, and the size of the droplet. These variables constitute two main dimensionless numbers, namely, the capillary number Ca = µc u/σ and the viscosity ratio λ = µd /µc , where u is the characteristic velocity of the flow, σ is the surface tension coefficient, and µd and µc are the viscosity of the dispersed and continuous phases, respectively. For very small sizes and slow velocities, the inertial effects could be neglected and so the Reynolds number (Re = ρc urd /µc where ρc is the density of the continuous phase and rd is the characteristic length scale) and the Weber 7 number (W e = ρc u2 rd /σ) are no longer important. The contribution of gravity is measured through the Bond number which is the ratio of the gravitational force to the surface tension force (Bo = ρ rd2 g/ σ). Generally micro-sized systems have a very small Bond number and therefore the Bond number is not an important dimensionless number. We will later see (by estimating the Bond number in chapter 3) that gravitational forces in our computations could be safely neglected. The first part of this study reports numerical simulations based on the Volume of Fluid method to investigate the influence of transmembrane pressure and crossflow velocity on the entry dynamics of thin oil films and droplets into pores of various cross-sections. We find that the formula derived in Ref. [65] for the critical pressure of permeation of an oil droplet into a circular pore agrees well with the results of numerical simulations. The numerical analysis is then extended to thin oil films covering pores with elliptical and rectangular cross-section in the absence of crossflow. In the presence of crossflow, we obtain numerically the phase diagram for the droplet rejection, permeation, and breakup as a function of the transmembrane pressure and shear rate, and study the details of the processes in three different regions of the phase diagram. These results are relevant to microfiltration of dilute oil-in-water emulsions at early stages before the formation of the cake layer. In the second part of this study, we perform numerical simulations to determine the drop behavior on a pore for different parameters. The effect of confinement on the droplet is demonstrated and an optimum height for the crossflow channel is found. Next, we perform simulations to study the effect of viscosity ratio, surface tension coefficient, contact angle, and drop-to-pore size ratio. Critical pressure of permeation is found and compared for different parameters. Also, breakup behavior of the droplet is studied. The rest of this report is organized as follows. The details of numerical simulations are 8 described in the next chapter. In chapter 3 analytical predictions based on the Young-Laplace equation are reviewed and verified numerically for an oil droplet on circular, elliptical, and rectangular pores. Moreover, the results for the oil droplet dynamics near circular pores in the presence of crossflow are presented. Chapter 4 introduces a novel procedure for computing the critical pressure of permeation and a summary of analytical predictions for the critical pressure based on the Young–Laplace equation is presented, and the effects of confinement, viscosity ratio, surface tension, contact angle, and droplet size on the critical transmembrane pressure and breakup are studied. The conclusions are given in the last chapter. 9 Chapter 2 Details of numerical simulations Numerical simulations were carried out using the commercial software FLUENT [94]. In order to control the transmembrane pressure and the crossflow velocity, a user-defined function was written in C and compiled along with the main solver. The Volume of Fluid method was used to solve the multiphase flow problem [95]. The Volume of Fluid method has the advantage that it rigorously conserves the mass and automatically captures topology changes. For a two-phase fluid, this method is based on the fact that the two phases form an impenetrable interface, i.e., each cell is filled with either one of the phases (denoting a specific phase zone) or a combination of two phases (denoting the interface). This is achieved by introducing a variable α, known as the “volume fraction”, which is defined as the ratio of the volume of fluid in the cell and the total cell volume; and it varies between 0 and 1 [96, 97]. An example of how the volume fraction varies near the interface is illustrated schematically in Fig. 2.1. The interface is tracked by solving the transport equation for the volume fraction as follows: ∂α + ∇⋅ (αV) = 0, ∂t (2.1) where V is the velocity vector. Equation (2.1) states that the substantial derivative of the volume fraction is zero, and, therefore, the interface is convected by the velocity fields at the interface. After solving Eq. (2.1), the material properties are computed by considering the 10 0.95 0.8 0.3 0 0 1 1 1 0.6 0 1 1 1 1 0.3 1 1 1 1 0.8 1 1 1 1 0.95 Figure 2.1 Values of the volume fraction for each of the phases and the interface (α = 1 in the first phase, α = 0 in the second phase, and 0 < α < 1 at the interface). fraction of each component in the cell; e.g., the density is given by ρ = α ρ2 + (1 − α) ρ1 , (2.2) where ρ is the volume-fraction-averaged density. One momentum equation is solved and the velocity field is shared between two phases as follows: ∂ (ρV) + ∇ ⋅ (ρVV) = −∇p + ∇ ⋅ [µ(∇V + ∇VT )] + ρ g + F, ∂t (2.3) where g is the vector of gravitational acceleration, and F is the source term. In multiphase flow applications, the source term is the surface tension force per unit volume and it is nonzero only at the interface. Using the divergence theorem, the surface tension force is defined 11 as the volume force in a cell as follows: ρ κ∇α Fσ = σ 1 , (ρ + ρ ) 1 2 2 (2.4) where σ is the surface tension between two phases and κ is the mean curvature of the interface in the cell. This force is related to the pressure jump across the interface (determined by the Young-Laplace equation) and it acts in the direction normal to the interface. The surface tension term tends to smooth out regions with large interface curvature [98]. If the interface is in contact with the wall, the normal vector (∇α), which defines the orientation of the interface in the cell adjacent to the wall, is determined by the contact angle. The effect of static contact angle is taken into account by imposing the interface unit normal for a point (a cell in the Finite Volume method), ni , on the wall containing the interface as follows: ni = nw cos(θst ) + nt sin(θst ), (2.5) where nw is the unit vector normal to the wall, nt is a vector on the wall and normal to the contact line, and θst is the static contact angle [76]. A SIMPLE algorithm was used for the pressure-velocity decoupling. The momentum equation was discretized using a second order upwind scheme. To reconstruct the interface and, consequently, solve the volume fraction transport equation, a PLIC (Piecewise Linear Interface Reconstruction) method was used [99]. The pressure equation is discretized using a staggered mesh with central differencing. In FLUENT, the interfacial tension is modeled using the well-known model of Continuum Surface Force (CSF) of Brackbill et al. [76]. Using CSF, the surface tension volume force Eq. (2.4) is added as a source term to the momentum 12 equation and the curvature is computed in terms of the vector normal to the interface n via: κ= 1 n [( ⋅ ∇)∣n∣ − (∇ ⋅ n)]. ∣n∣ ∣n∣ (2.6) Interfacial effects in multiphase flows are described by the Young-Laplace equation, which relates the pressure jump across the interface to its mean curvature and the surface tension coefficient. For flows at the micron length scales, the viscous effects are dominant and the inertial effects are typically negligible. The capillary number is a measure of how viscous shear stresses are compared to the interfacial stresses and it is defined Ca = µ U /σ, where U is the characteristic velocity, µ is the fluid viscosity, and σ is the surface tension coefficient. In chapter 3, the numerical simulations are performed to investigate a separation process of two commonly used liquids, i.e., kerosene and water. The density of kerosene is ρo = 889 kg/m3 and the viscosity ratio of kerosene and water at standard conditions is µo /µw = 2.4. It is assumed that the water is deionized; and, thus, the surface tension coefficient σ = 19.1 mN /m is used throughout the study [66]. Furthermore, we consider hydrophilic surfaces (for example polyvinyl-pyrrolidone [66]) with contact angles of kerosene in water greater than 90○ . In our simulations, the mesh was generated in GAMBIT using the Cooper mesh scheme. This method works by sweeping the node patterns of specified source faces through the whole volume and the resultant mesh consists of an array of hexagonal grids. For the results reported in the current study, we used about 30 cells along the pore diameter. To test the grid-resolution dependence, we considered 3 finer meshes that contained 50, 70, and 90 cells along the pore diameter. In the absence of crossflow, the simulations were performed for an oil droplet (rd = 1.0 µm) on a circular pore (rp = 0.2 µm) at two transmembrane pressures (1.000 and 0.951 bar) slightly above and below the exact value of the permeation 13 pressure 0.976 bar predicted by the Young-Laplace analysis. In all cases, the droplet would either penetrate into the pore or reside at the pore entrance for at least 40 µs. Furthermore, it was previously shown that the velocity of the contact line in the VOF method is inversely proportional to the logarithm of the mesh size [100], and, therefore, it is expected that the droplet velocity in the shear flow will depend on the grid resolution. However, in the present study, the oil droplet becomes temporarily pinned at the pore entrance by the transmembrane pressure, and thus the contact line velocity becomes much smaller than the flow velocity in the channel. Nevertheless, we have performed numerical simulations in the permeation, rejection, and breakup regions of the phase diagram with 4 times finer meshes and found that the predicted behavior of the the droplet remained unchanged. The corresponding profiles of the droplets simulated with the base mesh resolution level (30 grids along the diameter of the pore) and two and four times finer mesh are plotted in Fig. 2.2. There is great agreement between the profiles of two and four times refinement which indicates a grid-resolution-independent solution. To ensure that the numerical method used in this study correctly captures the physics of two-phase flows in stokes regime, we have performed a series of simulations of a well-studied benchmark problem in interface science, namely, the deformation and breakup of a droplet under pure shear flow. We have compared our results to findings from experiment, boundary integral method, and Volume of Fluid computations. The simulations were performed in a three-dimensional cubic space enclosed by periodic boundary conditions in lateral directions and no-slip walls in the vertical sides. The droplet is placed in the center of the cube and shear flow is induced by the motion of the walls in opposite directions. The viscosity and density of both fluids, and the shear rate is assumed to take a non-dimensional value of 1 and the drop radius is 0.25. Therefore, the Reynolds number is Re = ργr ˙ d2 /µ = 0.0625 which 14 -3E-06 X (µm) -2E-06 Base Level 2 Times Finer 4 Times Finer -1E-06 0 1E-06 -1E-06 0 1E-06 2E-06 3E-06 Y (µm) Figure 2.2 Profiles of a pinned droplet simulated using three levels of mesh resolution with the base level containing 30 grids along the diameter of the pore. allows for creeping flow. Fig. 2.3 shows the deformation and orientation angle of the droplet under various capillary numbers. It is observed that the present solver is able to accurately predict the deformation and the orientation angle of a single droplet in shear flow. The solver was able to numerically confirm that the critical capillary number for breakage of a droplet in simple shear flow is Cacr ≈ 0.41. We also comment that both the deformation and the orientation angle of the droplet appear to be increasing with the capillary number as will be further investigated in the following chapters. Moreover, Fig. 2.4 demonstrate the profiles of the droplet at four different capillary num- 15 Figure 2.3 Deformation (left) and orientation angle of a droplet in simple shear flow for various capillary numbers: ● present simulations; ◇ VOF simulations of Li et. al; ◻ boundary integral simulations of Rallison; × boundary integral simulations of Kwak et. al; △ boundary integral simulations of Kennedy et. al; ○ experimental results of Rumschiedt and Mason; the solid line is for small deformation theory of Cox [78] (parameters: ρ1 = ρ2 = 1, µ1 = µ2 = 1, rd = 0.25, and γ˙ = 1) 16 Ca 0.1 0.2 0.3 0.4 Figure 2.4 Snapshots of the cross-section of a droplet in simple shear flow for various capillary numbers using present simulations (left) and the results of VOF simulations of Li et. al (right) [78] (parameters: ρ1 = ρ2 = 1, µ1 = µ2 = 1, rd = 0.25, and γ˙ = 1) bers using the present solver and the numerical results of Li et. al [78]. There is remarkable resemblance in the shape of the droplet. The last shape corresponds to the capillary number slightly below the breakup value. 17 2 µm Y X Figure 2.5 Schematic of the hybrid mesh cross-section used in chapter 4. The mesh consists of coarse tetrahedral and fine hexagonal cells. In chapter 4, in order to increase the simulation efficiency, we generated a hybrid mesh that consists of fine hexagonal meshes in a part of the channel that contains the droplet and coarse tetrahedral meshes in the rest of the channel. Fig. 2.5 shows an example of the mesh used in chapter 4. Finally, a user-defined function was used to initialize the droplet shape and to adjust the velocity of the top wall that induced shear flow in the channel, as shown schematically in Fig. 2.6. Since the dynamics of the oil-water interface inside the pore slows down significantly when the transmembrane pressure becomes close to the critical pressure of permeation [101], the interface inside the pore is nearly static and the pressure jump across the spherical interface can be easily computed from the Young-Laplace equation. On the other hand, numerical simulations are required to resolve accurately the velocity fields, pressure, and the shape of the deformed droplet above the pore entrance. In the present study, we propose a new numerical procedure to compute the droplet critical pressure of permeation and breakup, as illustrated in Fig. 2.7. First, the pressure jump across the static interface inside the pore 18 Pressu re Inlet Periodic BC C Periodic B Moving Wall Symmetry Z Wall Y X Figure 2.6 Schematic representation of the oil droplet residing at the pore entrance in a rectangular channel with the corresponding boundary conditions. The width and length of the computational domain are fixed to 24 µm and 36 µm and the height of the domain is determined in Section 4.2. is calculated using the Young-Laplace equation. Second, we simulate the oil droplet in the presence of steady shear flow when the droplet covers the pore entrance completely and the oil phase fills in the pore. In the computational setup, the pore exit is closed to prevent the mass flux and to keep the droplet at the pore entrance. The difference in pressure across the deformed oil-water interface with respect to the inlet pressure is measured in the oil phase at the bottom of the pore (see Fig. 2.7). The critical pressure of permeation is then found by adding the pressure differences from the first and second steps. We applied this procedure to the cross-flow microflitration problem of chapter 3, in which the critical pressure for every shear flow is found by testing several operating conditions with similar shear rates and different transmembrane pressures. The present procedure was able to directly compute the critical pressures in a much shorter time and with a higher accuracy for various shear 19 P1 P1 = P2 P2 + P2 P3 P1 - P3 P3 = P1 - P2 + P2 - P3 Figure 2.7 Schematic of the droplet cross-section on the pore. Critical pressure of permeation (P1 − P3 ) is calculated in three steps: (1) pressure jump across static interface (P1 − P2 ) is calculated from Young-Laplace equation, (2) pressure jump across the dynamic interface (P2 − P3 ) is computed numerically, and (3) pressure jump from steps 1 and 2 are added. rates. Moreover, this procedure has the potential to be automated and requires less post processing. The solution of the Navier-Stokes equations for the flow over the membrane surface requires specification of the appropriate boundary conditions. As shown in Fig. 2.6, there are four types of boundary conditions used in the computational domain. The membrane and the pore surfaces are modeled as non-penetrable and no-slip boundaries. A moving “wall” boundary condition is applied at the top surface to induce shear flow between the moving top wall and the stationary membrane surface. The bottom of the pore is also described by the wall boundary condition to prevent mass flux and to keep the droplet pinned at the pore entrance (Note that, in the second part of chapter 3, the bottom of the pore is “pressure outlet” condition to allow flux and the imposition of transmembrane pressure. See Fig. 3.6). Periodic boundary conditions are imposed at the upstream and downstream of the channel. On the lateral side of the channel in the (Z+) direction, a pressure-inlet boundary condition 20 is applied to allow mass transfer, and to ensure that the reference pressure is fixed. Finally, a “symmetry” condition is implemented on the other side and only half of the computational domain is simulated to reduce computational efforts. We performed test simulations with an oil droplet rd = 2 µm exposed to shear flow and found that the local velocity profiles at the upstream, downstream, and the lateral sides remained linear when the width and length of the computational domain were fixed to 24 µm and 36 µm, respectively. The effect of confinement in the direction normal to the membrane surface on the droplet deformation and breakup will be investigated separately in the Section 4.2. 21 Chapter 3 Effect of Crossflow Velocity and Transmembrane Pressure on Microfiltration of Oil-in-Water Emulsions 3.1 The Young-Laplace analysis for circular pores Effective separation of oil-in-water emulsions is controlled by several key parameters such as the membrane pore size, surface energy, size of oil droplets, surface tension, and pressure difference across the membrane. If the transmembrane pressure is relatively high, then oil droplets will most likely penetrate the membrane resulting in low rejection rates of the oil phase. On the other hand, low transmembrane pressures tend to limit flux of water through the membrane. Hence, the optimum operating conditions strongly depend on the critical transmembrane pressure required for an oil droplet entry into a membrane pore. When the transmembrane pressure across a hydrophilic membrane exceeds a certain critical value, the oil phase will penetrate the membrane. Thus, for high separation efficiency, the transmembrane pressure should be maintained at a value below Pcrit , which for a continuous 22 oil film on the membrane surface with circular pores is given by Pcrit = 2 σ cos θ , rp (3.1) where σ is the surface tension coefficient between oil and water, θ is the contact angle of the interface of an oil droplet on a membrane surface immersed in water, and rp is the membrane pore radius [65]. The critical pressure in Eq. (3.1) is determined by the YoungLaplace pressure due to the curvature of the oil-water interface inside the pore. If instead of a thin oil film, a droplet of oil is placed at the entrance of the membrane pore, the formula for the critical pressure, Eq. (3.1), has to be corrected by a factor that depends on the ratio rd /rp to include the effect of the oil-water interface curvature above the membrane surface. It was previously shown [66, 65] that the pressure required to force an entry of an oil droplet of radius rd into a circular pore is given by Pcrit = 2 σ cos θ 3 2 + 3 cos θ − cos3 θ 1− . rp 4 (rd /rp )3 cos3 θ − (2 − 3 sinθ + sin3 θ) (3.2) In the limit rd → ∞, the curvature of the droplet above the pore vanishes, and thus this formula corresponds to a continuous oil film on the membrane surface, and Eq. (3.2) converges to Eq. (3.1). Contrary to the case of the thin film, the critical pressure for an oil droplet with θ = 90○ is negative, and the droplet will penetrate into the pore in the absence of the applied pressure gradient. We also comment that no such expression exists for the case when crossflow is present and the hydrodynamic drag force is exerted on the droplet parallel to the membrane surface. In what follows, we investigate the dynamics of an oil droplet and thin film entry into 23 pores with various cross-sections. The critical pressures [Eq. (3.1) and Eq. (3.2)] are compared with the results of numerical simulations in order to validate the numerical scheme. The numerical analysis is then extended to pores with rectangular and elliptical cross-sections. Finally, the effect of shear flow on the droplet entry is considered and a phase diagram of the transmembrane pressure versus shear rate is determined numerically. 3.2 An oil droplet on a circular pore in the absence of crossflow We first present the numerical results for the critical pressure required to force an entry of an oil droplet into a cylindrical pore. In the numerical scheme, the oil droplet is initially immersed in water above the membrane surface in the absence of flow. The transmembrane pressure is then gradually increased to a predetermined value. As the simulation continues, the droplet approaches the membrane surface and resides at the pore entrance. Depending on the applied pressure difference across the membrane, the droplet will either remain at the pore entrance or penetrate into the pore. We comment that when the applied pressure is close to the critical pressure, the dynamics of an oil droplet entry into the pore is significantly slowed down because the net driving force on the droplet is reduced. The critical pressure as a function of the droplet radius is plotted in Fig. 3.1 using Eq. (3.2) for three values of the pore radius. The error bars in Fig. 3.1 indicate the upper and lower values of the transmembrane pressure when the oil droplet either enters the pore or remains at the pore entrance during the time interval of about 40 µs. We find an excellent agreement between the results of numerical simulations and analytical predictions of Eq. (3.2), which provides a validation of the numerical method. Furthermore, as shown in Fig. 3.1, the critical 24 WATER 2 OIL Critical Pressure (bar) θ θ ∆P 1.5 × × × × 1 × × × × × × × × m rpore = 0 .2 µ × 0.5 × µm rpore = 0 .1 5 × × m rpore = 0 .3 µ × 0 0.4 0.6 0.8 1 1.2 1.4 Drop Radius (µm) Figure 3.1 The critical pressure computed using Eq. (3.2) for the values of the pore radius 0.15 µm (continuous line), 0.2 µm (dashed line), and 0.3 µm (dotted line). Error bars are extracted from the numerical simulations (see text for details) with the parameters µo /µw = 2.4, ρo /ρw = 0.781, σ = 19.1 mN /m, and θ = 135○ . The symbols indicate (×) rejection and (○) permeation of the oil droplet. pressure increases with increasing droplet radius or decreasing pore radius. As the droplet radius increases, the critical pressure approaches an asymptotic value predicted by Eq. (3.1) for a thin oil film (not shown). 25 3.3 Critical pressure for pores with arbitrary crosssection We next investigate the influence of pore cross-sectional shape on the critical transmembrane pressure using simple physical arguments and numerical simulations. According to the Young-Laplace equation, the pressure jump across an interface between two immiscible fluids is related to its mean curvature and the surface tension as follows: ∆P = 2 σ κ, (3.3) where κ is the mean curvature of the interface computed by averaging two principle curvatures. In the absence of gravity, the mean curvature of an arbitrary surface z(x, y) is given by 2 κ = ∇ ⋅ (√ ∇z 1 + ∣∇ z∣2 ). (3.4) It follows from Eq. (3.3) that an interface, which is subject to a prescribed pressure jump and constant surface tension coefficient, has a constant mean curvature. Therefore, if the gravity is negligible, the fluid-fluid interface forms a section of the so-called Delaunay surface [102, 103]. As shown in Fig. 3.2, if the oil-water interface is bounded by the walls of a pore of arbitrary cross-section, the constant contact angle at the pore surface imposes a boundary condition for Eq. (3.4) in the form cos θ = n ⋅ ( √ 26 ∇z 1 + ∣∇ z∣2 ), (3.5) C A OIL Interface z(x,y) n θ WATER y x Figure 3.2 Schematic representation of the oil-water interface inside the pore of arbitrary cross-section with perimeter C and area A. The interface forms a constant angle θ with the inner surface of the pore. where n is the outward unit vector normal to the pore surface [104]. Integrating Eq. (3.4) over an arbitrary cross-section with a smooth boundary and using the divergence theorem, we obtain Ap ⋅ 2 κ = ∫ Cp n ⋅ (√ ∇z 1 + ∣∇ z∣2 ) dL, (3.6) where Cp and Ap are the perimeter and cross-sectional area respectively. Taking the integral on the right hand side over the perimeter gives the following relation for the mean curvature of the oil-water interface [105] 2κ = Cp cos θ . Ap (3.7) According to Eq. (3.7), the mean curvature of an interface bounded by a pore of arbitrary 27 cross-section can be related to the geometric properties of the boundary. Remember that Eq. (3.3) relates the mean curvature to the surface tension and the pressure jump. In the case of an oil film on a pore with an arbitrary cross-section, the critical applied pressure is equal to the pressure jump at the interface. Therefore, combining Eqs. (3.3) and (3.7), we obtain the critical permeation pressure for the oil film to enter into a pore of arbitrary cross-section Pcrit = σ Cp cos θ . Ap (3.8) This equation can also be derived from the force balance between the applied pressure and the Laplace pressure due to the curvature of the oil-water interface inside the pore. It should be noted that the boundary value problem Eq. (3.4) does not always have a solution for a stable interface with a constant mean curvature and a constant contact angle [106]. In other words, there is a limitation on the values of the contact angle that correspond to the attached interface with a constant mean curvature. For example, if the contact angle (computed from the oil phase) is larger than the critical value, then Eq. (3.4) subject to the boundary condition Eq. (3.5) does not have a stable solution. As a result, the interface cannot remain attached to the bounding surface with a prescribed contact angle and, at the same time, maintain a constant mean curvature required by the Laplace equation. The critical contact angle is determined by the largest curvature of the crosssectional shape for smooth boundaries and by the smallest wedge angle for boundaries with sharp corners [107]. In what follows, we consider two special cases of rectangular and elliptical pores and compare predictions of Eq. (3.8) with the results of numerical simulations. The problem is illustrated schematically in the inset of Fig. 3.3. The oil film covers the pore entrance of a 28 hydrophilic membrane subject to a pressure gradient (the transmembrane pressure). Since the Bond number is small, Bo = (ρw − ρo ) rd2 g/ σ = 6 × 10−8 , the effect of the gravitational force can be neglected. 3.3.1 Thin oil film on the membrane surface with a rectangular pore In this subsection, we investigate the dynamics of an oil film entry into a rectangular pore, which is sometimes referred to as a “slotted pore” [108]. The rectangular shape of a slotted pore is characterized by the width (the shorter side) and the length (the longer side). Using Ap = w l and Cp = 2 (w + l), the corresponding critical pressure is obtained from Eq. (3.8) as follows: Pcrit = 2 σ cos θ ( 1 1 + ), w l (3.9) where w and l are the width and length of the pore cross-section [105, 109, 110]. In the limit when l ≫ w, Eq. (3.9) reduces to Pcrit = 2 σ cos θ / w, which is the critical pressure on an oil film entering into an infinitely long rectangular pore. In this case, one of the curvatures of the interface is zero and the other curvature is proportional to the width of the pore, and thus the shape of the interface is a part of a cylinder with the radius w/2 cos θ. The results of numerical simulations and predictions of Eq. (3.9) are shown in Fig. 3.3 for several aspect ratios. As expected, square pores have the highest critical pressure due to the largest perimeter-to-area ratio. The critical pressure decreases with increasing aspect ratio. These results demonstrate that there is an excellent agreement between the numerical results and analytical predictions based on the Young-Laplace equation. It is important to note that when the applied pressure is close to the critical pressure, the net force on the 29 WATER Critical Pressure (bar) 1 × 0.8 ∆P × × × × 0.6 × × 0.4 × l × × × × × × 1.2 1.4 × w = 0.4 µm × w = 0.6 µm × w = 1.0 µm w 0.2 0 OIL θ 0.4 0.6 0.8 1 1.6 l→ ∞ Length l (µm) Figure 3.3 The critical permeation pressure for the oil film into the rectangular pore with different aspect ratios. The curves are Eq. (3.9) and the symbols are the numerical results for µo /µw = 2.4, ρo /ρw = 0.781, σ = 19.1 mN /m, and θ = 120○ . The symbols indicate (×) rejection and (○) permeation of the oil film. interface is small, and, therefore, very long simulation time is required to capture the motion of the interface. The symbols in Fig. 3.3 indicate the rejection and permeation pressures that were resolved numerically without excessive computational effort. Interestingly, each curve in Fig. 3.3 is well described by the function Pcrit = 2 σ cos θ / l, which is shifted upward by a constant 2 σ cos θ / w (indicated by the horizontal lines in Fig. 3.3). It was previously shown that for a rectangular cross-section of arbitrary aspect ratio, Eq. (3.4) has a solution with a constant curvature for a non-wetting fluid (θ > 90○ ) when the the contact angle θ ≤ 135○ [111, 107, 112]. In the case of a square pore, the interface is part of a sphere with the radius w/2 cos θ [113, 104]. For other aspect ratios, the interface surface has a constant mean curvature κ = cos θ (1/w + 1/l), but it is no longer spherical [105]. 30 (a) θ = 90 o (b) θ = 120 o (c) θ = 135 o (d) θ = 150 o Figure 3.4 Snapshots of the oil-water interface inside the square pore for contact angles θ = 90○ , 120○ , 135○ , and 150○ . The critical pressure Eq. (3.9) is computed for the surface tension σ = 19.1 mN /m and the pore width 1.5 µm. Cases (a), (b), and (c) correspond to the stationary interface, while in the case (d) the interface is in transient state (see text for details). 31 Figure 3.4 shows snapshots of the oil-water interface inside the square pore obtained from our numerical simulations. The transmembrane pressure is set to a value computed from Eq. (3.9) for the contact angles θ = 90○ , 120○ , 135○ , and 150○ . As observed in Fig. 3.4, the concave shape of the interface strongly depends on the contact angle. When θ = 90○ , the interface enters the pore with zero curvature and, according to Eq. (3.9), with zero pressure gradient. For contact angles between 90○ and 135○ , the interface bends in the center and penetrates into the pore before it starts to move near the corners. For θ > 135○ , the distance between the interface location in the center and at the corners will theoretically be infinity because the interface becomes pinned at the corners while the inner part penetrates into the pore [114]. 3.3.2 Thin oil film on the membrane surface with an elliptical pore The critical permeation pressure for an oil film covering a pore with an elliptical cross-section can be estimated from Eq. (3.8) and the geometric properties of an ellipse. However, there is no exact expression for the perimeter of an ellipse. In our study, we use one of the most accurate and compact approximations that predicts the perimeter of an ellipse with an error of −0.04% [115] Cp ≈ π (a + b)[1 + 3h √ ], 10 + 4 − 3 h (3.10) where a and b are the major and minor radii of the ellipse and h = (a − b)2 /(a + b)2 . Using Eq. (3.10) and the expression for the ellipse area Ap = π a b, the critical pressure is given by Pcrit ≈ (a + b) 3h √ [1 + ] σ cos θ. ab 10 + 4 − 3 h 32 (3.11) WATER Critical Pressure (bar) 1 × ∆P × 0.8 × × × 0.6 OIL θ × × 0.4 × × × × × × × × × × 0.7 0.8 b= 0.2 µm b= 0.3 µm b= 0.5 µm 0.2 0 a 0.2 0.3 b 0.4 0.5 0.6 a→ ∞ Major Radius a (µm) Figure 3.5 The critical pressure of permeation for the oil film into the elliptical pore as a function of the major and minor radii. The curves are computed using Eq. (3.11). The symbols represent numerical results for the parameters µo /µw = 2.4, ρo /ρw = 0.781, σ = 19.1 mN /m, and θ = 120○ . The symbols indicate (×) rejection and (○) permeation of the oil film. Clearly, in the case of a circular pore, a = b = rp , Eq. (3.11) is reduced to Eq. (3.1). The results of numerical simulations and predictions of Eq. (3.11) are summarized in Fig. 3.5 for different aspect ratios. Similar to the case of the rectangular pore, the symbols indicate pressures of rejection and permeation that were resolved during the simulation time interval of about 40 µs. It can be observed that the critical pressure decreases with increasing ellipse aspect ratio, and the numerical results agree very well with predictions of Eq. (3.11). We also comment that when a ≫ b, each curve in Fig. 3.5 asymptotes to Pcrit ≈ 4 σ cos θ/π b, which is higher than the value Pcrit = 2 σ cos θ / w estimated for an infinitely long rectangular pore (see section 3.3.1). This difference arises because an infinitely long ellipse and an 33 infinitely long rectangle of equal width have the same perimeter but different areas. If the aspect ratio of an elliptical pore is less than 1.635, then the boundary value problem given by Eq. (3.4) has a solution for any contact angle [116]. However, when a/b > 1.635, there is a critical contact angle above which Eq. (3.4) has no solution [116]. For a/b = 1.635, the critical contact angle is 180○ , and it decreases with increasing aspect ratio. The largest aspect ratio considered in the present study is a/b = 4.0 for which the critical contact angle is 153.05○ [116]. Therefore, the contact angle of 120○ used in our simulations generated an interface with a constant mean curvature described by Eq. (3.7). We comment that the oil-water interface inside the elliptical pore is not spherical because the boundary condition Eq. (3.5) is not satisfied at the intersection of a sphere and a cylinder with an elliptical crosssection. Similar to the rectangular cross-section, we found that at the critical pressure given by Eq. (3.11), the contact line is pinned at the antipodal points of the highest curvature of the ellipse, (x, y) = (±a, 0), while the rest of the interface penetrates into the pore. 3.4 Sheared droplet on the membrane surface with a circular pore We next examine the combined effect of the transmembrane pressure and crossflow velocity on the entry dynamics of an oil droplet into a circular pore. The computational setup is illustrated schematically in Fig. 3.6. The shear flow is induced by translating the upper wall with a constant velocity. In our simulations, the effective shear rate γ˙ is defined as the ratio of the upper wall velocity to the channel height. The relevant dimensionless numbers, the capillary and Reynolds numbers, are estimated to be Ca = µw γr ˙ d /σ ≤ 0.03 and Re = ρw γr ˙ d2 /µw ≤ 0.5. The width of the channel is chosen to be about 8 times larger than the 34 H ≈ 2.3 rd Periodic BC Periodic BC Vw θ W L ≈ 10 rd rd 8 ≈ Pressure Outlet Figure 3.6 Schematic representation of the oil droplet in the channel with the circular pore. The shear flow is induced by the upper wall moving with a constant velocity parallel to the stationary lower wall. The droplet radius is rd = 0.9 µm and the pore radius is rp = 0.2 µm. Periodic boundary conditions are applied at the inlet and outlet, while a constant pressure is maintained at the side walls. droplet radius in order to minimize finite size effects in the lateral direction. Initially, the droplet is released upstream to insure that the flow reaches a steady state before the droplet approaches the pore. At the same time, the transmembrane pressure is set to a prescribed value, and the simulation continues until the droplet either reaches the outlet, penetrates into the pore, or breaks up. The main results of this study are summarized in Fig. 3.7, which shows the phase diagram for the droplet rejection, permeation, and breakup depending on the transmembrane pressure and shear rate. The corresponding snapshots of the droplet for five different cases (denoted by the capital letters A, B, C, D, and E) are presented in Fig. 3.8. Below, we discuss the details of the processes in the three different regions of the phase diagram and provide an estimate of the leakage volume during the droplet breakup. 35 Transmembrane Pressure (bar) E 1.15 Permeation 1.1 Breakup 1.05 B D A C 1 0.95 Rejection 0.9 0.85 1 2 3 4 5 6 Shear Rate (10 5 1/s) Figure 3.7 The phase diagram of the transmembrane pressure versus shear rate for the oil droplet with rd = 0.9 µm and the circular pore with rp = 0.2 µm. The contact angle is θ = 135○ . Each symbol represents a separate simulation that corresponds to either (◯) permeation, (square) breakup or (∇) rejection. Letters A, B, C, D, and E indicate operating conditions for the series of snapshots shown in Fig. 3.8. In the permeation region shown in Fig. 3.7, the transmembrane pressure is larger than the streamwise drag force, and the oil droplet penetrates into the pore. A series of snapshots at point B in Fig. 3.8 demonstrate the details of the permeation process. As observed in Fig. 3.8 (B), at first, the droplet partially penetrates into the pore and becomes strongly deformed in the shear flow. However, the droplet does not breakup because its size above the pore decreases as the droplet penetrates into the pore, and the viscous shear stress acts 36 t1 t2 t3 t4 t5 t6 t7 t8 A B C D E Figure 3.8 Sequences of snapshots of the sheared droplet on the pressurized pore for different operating conditions as indicated in Fig. 3.7. The droplet radius is rd = 0.9 µm, the pore radius is rp = 0.2 µm, and the contact angle is θ = 135○ . The letters denote (A) rejection at low shear rates (∆t ≈ 15 µs), (B) permeation (∆t ≈ 15 µs), (C) rejection at high shear rates (∆t ≈ 2 µs), (D) local breakup (∆t ≈ 2.5 µs), and (E) breakup with necking (∆t ≈ 4 µs). on a progressively smaller surface area. With increasing shear rate, the effect of viscous forces becomes more important, resulting in strong deformation of the droplet shape near the pore entrance. We find that at sufficiently large transmembrane pressures and Ca ≥ 0.015, the oil droplet breaks up. In this case, the larger droplet is washed off downstream and the smaller droplet enters the pore, and, as a result, the membrane leaks. Depending on the values of the transmembrane pressure and shear rate, two different breakup regimes were observed. The first regime is bounded by the minimum breakup pressure, which is found from Fig. 3.7 to be 1.00 ± 0.01 bar. Above this pressure, a small fragment is detached and penetrates into the pore while the main droplet is carried away by the shear flow. During this process, the droplet has a limited time to deform, and, therefore, the breakup occurs only locally without significant deformation, as shown in Fig. 3.8 (D). 37 The breakup process at higher transmembrane pressures occurs in a qualitatively different way [see Fig. 3.8 (E)]. As the droplet approaches the pore, it momentarily slows down and remains at the pore entrance during the “residence time”. In this case, the effects of the drag and the transmembrane pressure are relatively large and comparable with each other. As a result, the shape of the droplet is significantly deformed by the shear flow and a thin bridge is formed between two parts of the stretched droplet. This thinning is known as necking and it usually indicates the initial stage of the breakup process [117, 118]. For a short time interval, the thin neck holds the two parts of the droplet together. At the final stage of breakup, the neck gets thinner and thinner near the the edge of the pore, and at some point, it becomes unstable and the droplet breaks. Visual inspection of the snapshots of the droplet near the pore entrance revealed that, at a given shear rate in the breakup regime in Fig. 3.7, the residence time is roughly independent of the transmembrane pressure. Therefore, it is expected that during the breakup process, the volume of leaked droplets is proportional to the applied pressure. Figure 3.9 shows the leakage volume as a function of the transmembrane pressure when γ˙ = 5 × 105 s−1 . Indeed, the leakage volume is almost linearly proportional to the applied pressure, indicating that the flow inside the pore is described by the Hagen-Poiseuille equation. The lower part of the phase diagram in Fig. 3.7 indicates operating conditions when the oil droplet is rejected by the membrane and washed off by the shear flow. An example of the rejection process at low shear rates is presented in Fig. 3.8 (A). Although the droplet partially penetrates into the pore, the flow generates a force on the droplet surface, which pulls the droplet away from the pore, resulting in the droplet rejection. As shown in Fig. 3.8 (C), the residence time at higher shear rates is reduced, and the droplet is carried away by the flow without penetrating into the pore. 38 Transmembrane Pressure (bar) Rejected part 1.8 Leaked part 1.6 1.4 Numerical Simulation Best Linear Fit 1.2 1 0.05 0.1 0.15 0.2 3 Volume of Leakage (µm ) Figure 3.9 The leakage volume as a function of the applied pressure for rd = 0.9 µm, rp = 0.2 µm, and γ˙ = 5 × 105 s−1 . The square symbols indicate the numerical results and the dashed line is the best fit to the data. The inset shows a snapshot of the process shortly after breakup of the droplet. The contact angle of the oil droplet in water is θ = 135○ . The threshold of permeation is determined by the competition between the drag force and the transmembrane pressure. Naturally, with increasing shear rate, the drag force increases; and, therefore, it is not surprising that the boundary curve separating the permeation and rejection regions in Fig. 3.7 increases with shear rate. However, it is difficult to estimate the exact dependence of the drag force on the droplet because its shape becomes strongly deformed in the shear flow. We also comment that, in the range of shear rates reported in 39 Fig. 3.7, the lift force is about an order of magnitude smaller than the drag force [57]. The critical shear rate that marks the boundary between the permeation and breakup regions in Fig. 3.7 can be estimated using simple force balance arguments. In the absence of gravity, the sheared droplet is subject to forces of surface tension, Laplace pressure, drag, and lift. Neglecting the lift force and following the analysis in Ref. [57], the torque balance around the edge of the pore for the droplet configuration depicted in Fig. 3.8 (E t6 ) can be written as follows: FD dd + (P2 − P1 ) An dn − Fσ dn = 0, (3.12) where dn and An = π dn2 /4 are the diameter and cross-sectional area of the thinnest stable neck, Fσ = π σ dn is the surface tension force around the perimeter of the neck, and P2 and P1 are the pressures inside the droplet and in the channel, respectively. The Stokes drag force on the spherical droplet, FD ≈ 1.1 π µw γd ˙ d2 , is estimated for the viscosity ratio 2.4 and the average flow velocity γd ˙ d /2 [58]. In our simulations, the typical diameter of the thinnest stable neck is dn ≈ 0.9 dp [see Fig. 3.8 (E t6 )]. Using dd = 1.7 µm, σ = 19.1 mN /m, and µw = 10−3 kg/m s in Eq. (3.12), the critical shear rate is roughly estimated to be γ˙ ≈ 3.4 × 105 s−1 , which is in good agreement with the value γ˙ ≈ 3.2 × 105 s−1 obtained numerically in Fig. 3.7. The permeation, rejection, and breakup regions identified in the phase diagram in Fig. 3.7 can be useful for the optimal design and operation of crossflow microfiltration systems. It is apparent that the permeation region should be avoided for filtration purposes. The optimal performance of the microfiltration system with maximum rejection is achieved in the upper part of the rejection region where the large transmembrane pressure results in high flux of water while the oil phase is completely rejected. However, the separation efficiency can be increased at higher transmembrane pressures in the breakup region, where the higher flux 40 of water is accompanied by some oil leakage. 41 Chapter 4 Effect of Physical and Geometrical Parameters on Crossflow Microfiltration of Liquid-Liquid Emulsions 4.1 Analytical formulations As discussed in Section 3.3, The pressure jump across a static interface between two immiscible fluids could be described in terms of the interfacial tension coefficient σ and the mean curvature of the interface κ by the Young-Laplace equation as ∆P = 2 σ κ. Therefore, the pressure required to force entry (i.e. critical pressure of permeation) of a thin liquid film into a pore of arbitrary cross-section is found by Eq. 3.8 and was used in chapter 3 to calculate critical pressure of permeation of thin oil films into rectangular, elliptical and circular pores and compared to results obtained numerically. In case of a droplet on a pore, we need to account for the interface of the droplet inside the channel as well as the interface in the pore and it was shown that Eq. 3.2 could be used to find the critical pressure of permeation of a drop on a circular pore. 42 Numerical simulations of Section 3.2 for a drop size of 0.15 µm and different pore sizes show very good agreement with Eq. 3.2. In the presence of cross-flow in the channel, Eq. 3.2 no longer holds since the drop interface is distorted and assumption of a perfectly spherical interface is not valid. Furthermore, numerical simulations showed that the critical pressure of permeation increases with increasing crossflow velocity up to a certain value, above which the droplet breaks up. [101]. The region of critical pressure is denoted by the area between the permeation region and the rejection region in the phase diagram Fig. 3.7. In this chapter we calculate the critical pressure region more accurately in the form of a line and determine its dependence on shear rate for a range of material properties and geometrical parameters. In the presence of crossflow above the membrane surface, an oil droplet breaks up when viscous stresses over the droplet surface exposed to the flow become larger than capillary stresses at the interface of the droplet near the membrane pore. Therefore, at the moment of breakup, the drag force in the flow direction is balanced by the capillary force at the droplet interface around the pore D ≈ Fσ . (4.1) Neglecting the contact angle dependence, Fσ ∝ σ rp is the interfacial force acting in the direction opposite to the flow at the droplet interface near the pore entrance. The drag force generated by a linear shear flow on a spherical droplet attached to a solid surface is given by D ∝ fD (λ) µ γ˙ rd2 , (4.2) where µ is the viscosity of the continuous phase, γ˙ is the shear rate, and rd is the radius of the droplet [119, 85]. The coefficient fD (λ) is a function of the viscosity ratio λ = µoil /µwater and it depends on the shape of the droplet above the surface. Sugiyama and Sbragaglia [85] 43 have estimated this function analytically for a hemispherical droplet (θ = 90○ ) attached to a solid surface fD (λ) ≈ 2 + 4.510 λ . 1 + 1.048 λ (4.3) By plugging Eq. (4.2) into Eq. (4.1) and introducing r¯ = rd / rp , the critical capillary number for breakup of a droplet on a pore can be expressed as follows: Cacr ∝ 1 , fD (λ) r¯ (4.4) where the capillary number is defined as Ca = µw γr ˙ d /σ. The difference in pressure inside the pore in the presence of flow and at zero shear rate can be estimated from the torque generated by the shear flow on the droplet surface. The torque around the center of the droplet projected on the membrane surface is given by T ∝ fT (λ) µ γ˙ rd3 , (4.5) It was previously shown [85] that for a hemispherical droplet on a solid surface, fT (λ) is a function of the viscosity ratio fT (λ) ≈ 2.188 λ . 1 + 0.896 λ (4.6) Hence, the balance of the torque due to shear flow above the membrane surface [given by Eq. (4.5)] and the torque arising from the pressure difference, (Pcr − Pcr0 ) Ap rd , can be reformulated in terms of the capillary number and drop-to-pore size ratio as follows: f (λ) σ r¯ Ca Pcr − Pcr0 ∝ T , rp 44 (4.7) where Pcr0 is the critical permeation pressure in the absence of crossflow. In what follows, we consider the effects of confinement, viscosity ratio, surface tension, contact angle, and droplet size on the critical pressure of permeation and breakup using numerical simulations and analytical predictions of Eq. (4.4) and Eq. (4.7). 4.2 The effect of confinement on the critical pressure of permeation and breakup In practical applications, the dimensions of a crossflow channel of a microfiltration system are much larger than the typical size of emulsion droplets so that the velocity profile over the distance of about rd from the membrane surface can be approximated as linear. To more closely simulate this condition in our computational setup, the shear flow above the membrane surface was induced by moving the upper wall of the crossflow channel (Fig. 2.6). To understand how the finite size of the channel affects droplet dynamics at the membrane surface, we studied the influence of the channel height on the droplet behavior. The confinement ratio is defined as the ratio of the height of the droplet residing on the pore at zero shear rate Hd (i.e., the height of a spherical cap above the membrane surface) to the channel height Hch . It is important to note that the degree of confinement is varied only in the direction normal to the membrane surface and the computational domain is chosen to be wide enough for the lateral confinement effects to be negligible (see chapter 2). We performed numerical simulations of an oil droplet with radius rd = 2 µm in steadystate shear flow for the channel heights 3.8 µm ≤ Hch ≤ 12.0 µm. Figure 4.1 illustrates the effect of confinement on the shape of the droplet residing on a rp = 0.5 µm pore when the capillary number is Ca = µw γr ˙ d /σ = 0.021. The height of the droplet above the membrane 45 Flow Direction 1µm Hd/Hch = 0.286 Hd/Hch = 0.428 Hd/Hch = 0.686 Hd/Hch = 0.902 Figure 4.1 The cross-sectional profiles of oil droplets in steady shear flow for the indicated confinement ratios when the capillary number is Ca = 0.021. The droplet radius is rd = 2 µm, the pore radius is rp = 0.5 µm, the contact angle is θ = 135○ , the surface tension coefficient is σ = 19.1 mN /m, and the viscosity ratio is λ = 1. surface in the absence of flow is approximately 3.43 µm. It can be observed from Fig. 4.1 that highly confined droplets become more elongated in the direction of flow than droplets with lower confinement ratios, which is in agreement with the results of previous simulations [120]. When a droplet is highly confined, the distance between the upper moving wall and the top of the droplet is relatively small. As a result, the effective shear rate at the surface of the droplet is higher and the droplet undergoes larger deformation. Furthermore, the crosssectional profiles for the confinement ratios of 0.428 and 0.286 are nearly identical, indicating that the flow around the droplet is not affected by the upper wall when Hd /Hch ≤ 0.428 and 46 0.032 Cacr 0.029 0.026 0.023 0.2 0.4 0.6 0.8 1 Confinement Ratio (Hd/Hch) Figure 4.2 The critical (breakup) capillary number as a function of the confinement ratio Hd /Hch . Other parameters are the same as in Fig. 4.1. the capillary number is fixed. Figure 4.2 shows the variation of the critical capillary number (right before breakup) as a function of the confinement ratio for the the same material parameters as in Fig. 4.1. These results indicate that highly confined droplets breakup at lower capillary numbers, and, when the confinement ratio is smaller than about 0.5, the breakup capillary number remains nearly constant. For the rest of the study, the channel height was fixed to 8 µm, which corresponds to the confinement ratio of 0.428 for a droplet with radius rd = 2 µm. For the results presented in the subsection 4.6, the channel height was scaled appropriately to retain the same confinement ratio for larger droplets. 4.3 The effect of viscosity ratio on the critical transmembrane pressure The ratio of viscosities of the dispersed and continuous phases is an important factor that determines the magnitude of viscous stresses at the interface between the two phases. For a 47 small droplet at low Reynolds numbers, the viscous stresses are primarily counterbalanced by interfacial tension stresses. In a shear flow, viscous stresses tend to distort the surface of a droplet, while interfacial stresses assist in retaining its initial spherical shape. The competition between the two stresses determines the breakup criterion, deformation, and orientation of the droplet [62, 121]. In this subsection, we investigate numerically the effect of viscosity ratio on the droplet deformation and breakup at the entrance of the membrane pore. Figure 4.3 shows the effect of the viscosity ratio, λ = µo /µw , on the critical pressure of permeation and breakup of an oil droplet on a membrane pore as a function of the capillary number. The percent increase in critical pressure is defined with respect to the critical pressure in the absence of crossflow Pcr0 , i.e., (Pcr − Pcr0 )/Pcr0 × 100 %. Keeping in mind that Pcr0 does not depend on λ, the results shown in Fig. 4.3 demonstrate that at a fixed Ca, the critical pressure increases with increasing viscosity ratio, which implies that higher viscosity droplets penetrate into the pore at higher transmembrane pressures. Specifically, the maximum increase in critical pressure just before breakup is about 8 % for λ = 1 and about 15 % for λ = 20. Furthermore, highly viscous droplets tend to break at lower shear rates because of the larger drag force generated by the shear flow [see Eq. (4.2)]. As reported in Fig. 4.3, the critical capillary number for breakup varies from about 0.018 for λ = 20 to 0.032 for λ = 1. The practical implication of these results is that in membrane emulsification processes the use of liquids with lower viscosity ratios should be avoided as the droplets tend to break at higher shear rates. Examples of cross-sectional profiles of the oil droplet in steady shear flow are presented in Fig. 4.4 for the viscosity ratio λ = 1. At small capillary numbers, no significant deformation occurs and the droplet retains its spherical shape above the membrane surface. As Ca 48 % Increase in Critical Pressure 16 12 8 λ=1 λ=2 λ=3 λ=5 λ = 10 λ = 20 4 0 0 0.007 0.014 0.021 0.028 0.035 Ca Figure 4.3 The percent increase in critical pressure of permeation as a function of the capillary number Ca = µw γr ˙ d /σ for the indicated viscosity ratios λ = µo /µw . Typical error bars are shown on selected data points. For each value of λ, the data are reported up to the critical capillary number above which droplets break into two segments. The droplet and pore radii are rd = 2 µm and rp = 0.5 µm, respectively. The contact angle is θ = 135○ and the surface tension coefficient is σ = 19.1 mN /m. increases, a neck forms at the pore entrance while the rest of the droplet remains nearly spherical. A closer look at the shapes of the droplet for Ca = 0.0283 and 0.0314 in Fig. 4.4 reveals that with increasing shear flow, the neck gets thinner and the droplet becomes more elongated in the direction of flow. While the torque due to the shear flow does not increase significantly, the elongated shape of the droplet results in an effectively longer arm for the 49 Flow Direction 1µm Ca = 0.0063 Ca = 0.0126 Ca = 0.0188 Ca = 0.0251 Ca = 0.0283 Ca = 0.0314 Figure 4.4 The cross-sectional profiles of the oil droplet residing on the circular pore with rp = 0.5 µm for the indicated capillary numbers. The viscosity ratio is λ = 1. Other parameters are the same as in Fig. 4.3. torque due to pressure in the droplet along the flow direction, and, thus, it leads to a lower critical permeation pressure required to keep the droplet attached to the pore. This effect is observed in Fig. 4.3 as the critical pressure just before breakup decreases as a function of Ca. We next estimate the breakup time and compare it with the typical deformation time of the droplet interface for different viscosity ratios. In our simulations, the upper wall velocity is increased quasi-steadily and the spontaneous initiation of the breakup process can be clearly detected by visual inspection of the droplet interface near the pore entrance. We then identify the moment when a droplet breaks into two segments and compute the breakup time. The deformation time scale, defined by µw rd (1 + λ)/σ, is a measure of the typical relaxation 50 200 Flow Direction Breakup Time (µs) 1µm 150 100 λ=1 λ=2 λ=3 λ=5 λ = 10 λ = 20 50 0 0 0.6 1.2 1.8 2.4 Deformation Time Scale, µwrd(1 + λ)/σ (µs) Figure 4.5 The breakup time versus deformation time scale µw rd (1 + λ)/σ for the tabulated values of the viscosity ratio λ = µo /µw . Other system parameters are the same as in Fig. 4.3. The straight line is the best fit to the data. The error bars for the breakup time are about the symbol size. The inset shows the droplet profiles just before breakup for the same viscosity ratios. time of the droplet interface with respect to its deformation at steady state [80, 122]. In Fig. 4.5, the breakup time is plotted against the deformation time scale for different viscosity ratios. Notice that the breakup time increases linearly with the deformation time scale, which confirms that highly viscous droplets break up more slowly. The inset in Fig. 4.5 displays the droplet cross-sectional profiles just before breakup for the same viscosity ratios. It can be observed that the profiles nearly overlap with each other, indicating that droplets 51 (% Increase in Critical Pressure) / fT(λ) 8 λ=1 λ=2 λ=3 λ=5 λ = 10 λ = 20 6 4 2 0 0 0.02 0.04 0.06 0.08 0.1 Ca fD(λ) Figure 4.6 The normalized percent increase in critical pressure of permeation versus the modified capillary number Ca fD (λ) for the selected values of the viscosity ratio λ = µo /µw . The functions fD (λ) and fT (λ) are given by Eq. (4.3) and Eq. (4.6), respectively. with different viscosities are deformed identically just before breakup. According to Eq. (4.4), the breakup capillary number depends on the drop-to-pore size ratio and the viscosity ratio via the function fD (λ). Therefore, it is expected that the product Cacr fD (λ) will be independent of λ and the appropriate dimensionless number for a constrained viscous droplet in a shear flow is Ca fD (λ). Moreover, based on Eq. (4.7), the percent increase in the critical pressure is independent of the viscosity ratio when it is divided by fT (λ). Figure 4.6 shows the same data as in Fig. 4.3 but replotted in terms of the 52 normalized critical pressure and the modified capillary number. As is evident from Fig. 4.6, the data for different viscosity ratios nearly collapse on the master curve. It is seen that droplets break at approximately the same value Ca fD (λ) ≈ 0.09. In practice, the increase in critical pressure due to crossflow can be roughly estimated from the master curve in Fig. 4.6 for any viscosity ratio in the range 1 ≤ λ ≤ 20. Also, if Ca fD (λ) ≥ 0.09, the oil droplets will break near the pore entrance for any viscosity ratio. 4.4 The effect of surface tension on the critical pressure of permeation Surface tension coefficient determines how well an interface can resist external forces such as pressure and shear stresses. In the absence of flow, the normal component of the surface tension force is canceled by a pressure difference across the interface known as capillary pressure which relates to the surface tension coefficient through the Young-Laplace equation. In the presence of shear stresses, however, the shape of the interface is also a function of the viscosity ratio and flow velocity. Here, we investigate the influence of surface tension on the critical permeation pressure, deformation and breakup of an oil droplet residing at the pore entrance in the presence of crossflow above the membrane surface. Figure 4.7 shows the critical pressure of permeation as a function of shear rate for five values of the surface tension coefficient. As expected from Eq. (3.2), the critical pressure at zero shear rate increases linearly with increasing surface tension coefficient. Note that oil droplets with higher surface tension break up at higher shear rates because larger stresses are required to deform the interface and cause breakup of the neck. Also, the difference between the critical pressure just before breakup and Pcr0 is larger 53 Critical Pressure of Permeation (kPa) 77 σ = 38.2 mN /m × 64 51× 38 × σ = 28.6 mN /m σ = 19.1 mN /m σ = 14.3 mN /m 25× σ = 9.55 mN /m × 12 0 1.2 2.4 3.6 4.8 6 5 Shear Rate (10 1/s) Figure 4.7 The critical pressure of permeation as a function of shear rate for the indicated surface tension coefficients. The symbols (×) denote the analytical predictions of Eq. (3.2). The droplet and pore radii are rd = 2 µm and rp = 0.5 µm, respectively. The viscosity ratio is λ = 1 and the contact angle is θ = 135○ . at a higher surface tension; for example, it is about 1.5 kPa for σ = 9.55 mN /m and 6 kPa for σ = 38.2 mN /m. The results shown in Fig. 4.7 suggest that crossflow microfiltration of emulsion droplets with higher surface tension is more efficient because higher transmembrane pressure can be applied and the droplet breakup is less likely. Examples of droplet cross-sectional profiles above the membrane pore are presented in Fig. 4.8 for five values of the surface tension coefficient. These profiles are extracted from 54 Flow Direction 1µm σ = 0.0382 N/m σ = 0.0286 N/m σ = 0.0191 N/m σ = 0.0143 N/m σ = 0.00955 N/m Figure 4.8 The cross-sectional profiles of the oil droplet above the circular pore for the listed values of the surface tension coefficient. In all cases, the shear rate is γ˙ = 1.5 × 105 s−1 . Other parameters are the same as in Fig. 4.7. the data reported in Fig. 4.7 at the shear rate γ˙ = 1.5 × 105 s−1 . It can be observed that oil droplets with lower surface tension become highly deformed along the flow direction. The elongation is especially pronounced when the surface tension coefficient is small; for σ = 9.55 mN /m the droplet interface is deformed locally near the pore entrance and the neck is formed. To further investigate the effect of surface tension on the droplet breakup, we compare the breakup time and the deformation time scale µw rd (1 + λ)/σ. The numerical results are summarized in Fig. 4.9 for the same values of the surface tension coefficient as in Fig. 4.7. Similar to the analysis in the previous subsection, the breakup time was estimated from the 55 35 Flow Direction Breakup time (µs) 28 1µm 21 14 σ = 38.2 mN/m σ = 28.6 mN/m σ = 19.1 mN/m σ = 14.3 mN/m σ = 9.55 mN/m 7 0 0 0.1 0.2 0.3 0.4 Deformation Time Scale, µwrd(1 + λ)/σ (µs) Figure 4.9 The breakup time versus deformation time scale µw rd (1 + λ)/σ for the surface tension coefficients in the range from 9.55 mN /m to 38.2 mN /m. Other parameters are the same as in Fig. 4.7. The straight solid line is the best fit to the data. The cross-sectional profiles of the oil droplet just before breakup are displayed in the inset. time when a droplet becomes unstable under quasi-steady perturbation till the formation of two separate segments. It can be observed in Fig. 4.9 that the breakup time varies linearly with increasing deformation time scale, which in turn indicates that the breakup time is inversely proportional to the surface tension coefficient. In addition, the inset in Fig. 4.9 shows the cross-sectional profiles of the droplet just before breakup for the same surface tension coefficients. Interestingly, the profiles nearly coincide with each other, indicating 56 % Increase in Critical Pressure 10 σ = 9.55 mN/m σ = 14.3 mN/m σ = 19.1 mN/m σ = 28.6 mN/m σ = 38.2 mN/m 8 6 4 2 0 0 0.007 0.014 0.021 0.028 0.035 Ca Figure 4.10 The percent increase in critical pressure of permeation as a function of the capillary number Ca = µw γr ˙ d /σ for the selected values of the surface tension coefficient. The rest of the material parameters are the same as in Fig. 4.7. that the droplet shape at the moment of breakup is the same for any surface tension. In order to present our results in a more general form, we replotted the data from Fig. 4.7 in terms of the percent increase in critical pressure, (Pcr − Pcr0 )/Pcr0 × 100 %, and the capillary number in Fig. 4.10. Note that in all cases, the data collapse onto a master curve and breakup occurs at the same relative pressure (Pcr − Pcr0 )/Pcr0 ≈ 8 % and Cacr ≈ 0.03, which indicates that the capillary number is an appropriate dimensionless number to describe the droplet deformation in shear flow with variable surface tension. These results are not 57 surprising, given that the breakup capillary number, Eq. (4.4), does not depend on the surface tension coefficient. Moreover, the increase in critical pressure due to crossflow, Eq. (4.7), is proportional to σ and Ca, and when it is divided by Pcr0 , which itself is a linear function of σ [see Eq. (3.2)], the percent increase in critical pressure becomes proportional to the capillary number. In practice, the master curve reported in Fig. 4.10 can be used to predict the critical permeation pressure and breakup of emulsion droplets for specific operating conditions and surface tension. 4.5 The effect of contact angle on permeation of the droplet The problem of a drop on a solid surface involves three phases (the drop, the continuous phase, and the solid surface) that intersect at a contact line with a specific contact angle. The contact angle is determined by the contribution of the surface tension forces between either on of the two phases and is a measure of the direction of the surface tension force on the at the contact line. Here, we focus on the effect of contact angle on the permeation pressure, deformation and breakup of oil droplets on a membrane pore. The variation of the critical permeation pressure as a function of the capillary number is presented in Fig. 4.11 for nonwetting oil droplets with contact angles 115○ ≤ θ ≤ 155○ . The critical pressure at zero shear rate is higher for oil droplets with larger contact angles, which is in agreement with the analytical prediction of Eq. (3.2). As expected, with increasing shear rate, the critical pressure of permeation increases for all values of θ studied. We estimate the maximum change in the critical pressure to be about 3 kPa and roughly independent of the contact angle. This 58 Critical Pressure of Permeation (kPa) 57 × 48 × 39 × 30 × θ = 155 o θ = 145 o θ = 135 o θ = 125 o 21 o θ = 115 × 12 0 0.6 1.2 1.8 2.4 3 5 Shear Rate (10 1/s) Figure 4.11 The critical pressure of permeation as a function of the capillary number for the indicated contact angles. The critical pressure at zero shear rate, given by Eq. (3.2), is denoted by the symbols (×). The droplet radius, pore radius, surface tension coefficient, and viscosity ratio are rd = 2 µm, rp = 0.5 µm, σ = 19.1 mN /m, and λ = 1, respectively. corresponds to a relative increase of about 6% for the contact angle θ = 155○ and 21% for θ = 115○ . These results suggest that the relative efficiency of a microfiltration system due to crossflow is higher for emulsion droplets with lower contact angles. Interestingly, we find that the critical capillary number for breakup (Cacr ≈ 0.032) is nearly independent of the contact angle. This suggests that Ca can be used as a criterion for predicting breakup. Finally, the examples of the droplet cross-sectional profiles are shown in Fig. 4.12 for different contact 59 Flow Direction 1µm θ = 115o θ = 125o θ = 135o o θ = 145 θ = 155o Figure 4.12 The cross-sectional profiles of the oil droplet above the circular pore for the listed values of the contact angle when Ca = 0.022. Other parameters are the same as in Fig. 4.11. angles when Ca = 0.022. Notice that droplets with lower contact angles wet larger solid area and are less tilted in the direction of flow. 4.6 The effect of droplet size on the critical pressure of permeation In the microfiltration process, the size of the membrane pore is one of the crucial parameters that determine the permeate flux and membrane selectivity. Membranes with smaller pore sizes provide higher rejections but require higher transmembrane pressures to achieve the 60 Critical Pressure of Permeation (kPa) 42 = r d/r p 39× × 36 5 4 = r / p .5 rd = rd/rp × 4 =3 rd/rp .5 × 33 rd/rp = 30 3 × 27 0 1.1 2.2 3.3 4.4 5.5 5 Shear Rate (10 1/s) Figure 4.13 The critical pressure of permeation as a function of shear rate for the selected drop-to-pore size ratios. The symbols (×) indicate the critical pressure in the absence of flow calculated from Eq. (3.2). The pore radius, surface tension coefficient, contact angle, and viscosity ratio are rp = 0.5 µm, σ = 19.1 mN /m, θ = 135○ and λ = 1, respectively. same permeate flux. Moreover, size of the droplets generated in a membrane emulsification system greatly depends on the size of the pores [40, 123]. Emulsions with narrower size distribution tend be more stable [124]. Therefore, it is essential for the droplet size and the pore size to match appropriately. In this subsection, we examine the influence of the drop-to-pore size ratio on the critical pressure of permeation and the breakup dynamics of oil droplets in the presence of crossflow above the membrane surface. 61 Figure 4.13 reports the critical permeation pressure as a function of shear rate for the droplet radii in the range from 1.5 µm to 2.5 µm, while the pore radius is fixed at rp = 0.5 µm. In the absence of crossflow, the critical pressure is higher for larger droplets because they have lower curvature of the interface above the membrane surface, which is in agreement with the analytical prediction of Eq. (3.2). With increasing shear rate, the critical pressure increases for droplets of all sizes. Note also that the slope of the curves in Fig. 4.13 is steeper for larger droplets because of the larger surface area exposed to shear flow, resulting in a higher drag torque, and, consecutively, a higher transmembrane pressure needed to balance the torque. Furthermore, as shown in Fig. 4.13, smaller droplets break at higher shear rates, since higher shear stress are required to produce sufficient deformation for the breakup to occur. The maximum relative critical pressure is about 14% for rd /rp = 3 and 6% for rd /rp = 5. We next compute the difference in the critical permeation pressure with respect to the critical pressure in the absence of flow, Pcr − Pcr0 , and define r¯ = rd /rp . According to Eq. (4.4), the product Cacr × r¯ is independent of the droplet radius. At the same time, Eq. (4.7) suggests that the increase in critical pressure depends on the droplet radius via the term Ca × r¯. Figure 4.14 shows the critical pressure difference as a function of the modified capillary number Ca × r¯ for different droplet radii. It can be observed in Fig. 4.14 that all curves nearly collapse on each other and the droplet breakup occurs at the same value Ca × r¯ ≈ 0.125. We also comment that one of the assumptions in deriving Eq. (4.7) is that the distance between the center of the pore and the center of the droplet on the membrane surface is approximately rd . This approximation becomes more accurate for larger drop-topore size ratios, and, thus, the critical pressure difference in Fig. 4.14 is nearly the same for larger droplets even at high shear rates. The inset of Fig. 4.14 shows the cross-sectional profiles of oil droplets just before breakup 62 5 Flow Direction 1µm 3 0 Pcr - Pcr (kPa) 4 2 rd/rp = 3.0 rd/rp = 3.5 rd/rp = 4.0 rd/rp = 4.5 rd/rp = 5.0 1 0 0 0.03 0.06 0.09 0.12 0.15 Ca ×r Figure 4.14 The difference in the critical pressure, Pcr − Pcr0 , versus the modified capillary number for five drop-to-pore size ratios r = rd /rp and rp = 0.5 µm. Other parameters are the same as in Fig. 4.13. The cross-sectional profiles of the droplet just before breakup are shown in the inset for same r. for different drop-to-pore size ratios. Note that all droplets are pinned at the pore entrance and elongated in the direction of flow. It is seen that when r¯ is small, the droplet shape is significantly deformed from its original spherical shape. In contrast, larger droplets remain nearly spherical and only deform near the pore entrance. In general, the droplet-to-pore size ratio should be large enough to make Pcr sufficiently high for practicable separation. At the same time, if the pore size is much smaller than the droplet size, the water flux through the 63 membrane decreases and the probability of breakup increases, which could result in lower rejection rates and internal fouling of the membrane. Therefore, choosing a membrane with an appropriate pore size could greatly increase the efficiency of the microfiltration process. 64 Chapter 5 Conclusions In this report, we investigated numerically the effects of material parameters and operating conditions on the entry dynamics of thin oil films and droplets into pores of a microfiltration membrane. Specifically, the critical pressure of permeation, breakup criteria, and deformation of oil droplets were studied for a number of material parameters, such as viscosity ratio, surface tension, contact angle, and drop size, and reported for different crossflow velocities above the membrane surface. We used a finite-volume-based incompressible flow solver that used Volume of Fluids method to track the interface between water and oil. The validity of the solver was confirmed by testing benchmark problems of droplets in simple shear flow. In the first part, the numerical simulations were carried out to investigate the influence of the transmembrane pressure and crossflow velocity on the entry dynamics of thin oil films and oil droplets into pores of different cross-section. We considered hydrophilic membrane surfaces with contact angles of oil in water greater than 90○ . The numerical method was validated against the analytical solution for the critical pressure of permeation of an oil droplet into a circular pore in the absence of crossflow. Furthermore, we found that the results of numerical simulations of thin oil films on elliptical or rectangular pores agree well with the theoretical prediction for the critical pressure expressed in terms of geometric parameters of the pore cross-section. Also, examples of curved oil-water interfaces inside elliptical and rectangular pores were discussed for different aspect ratios and contact angles. In the presence of crossflow above the membrane surface, we have determined numerically 65 the phase diagram for the droplet rejection, permeation, and breakup as a function of the transmembrane pressure and shear rate. A detailed analysis of the droplet dynamics near the pore entrance was performed in the three different regions of the phase diagram. We found that in the permeation region, the transmembrane pressure is larger than the streamwise drag, and the oil droplet penetrates into the circular pore. With increasing crossflow velocity, the shape of the droplet becomes strongly deformed near the pore entrance; and, at sufficiently high transmembrane pressures and shear rates, the droplet breaks up into two fragments, one of which penetrates into the pore. It was also shown that during the breakup process, the residence time of the oil droplet at the pore entrance is roughly independent of the transmembrane pressure, and the volume of the leaked fragment is nearly proportional to the applied pressure. Finally, the numerical value of the critical shear rate that separates the permeation and breakup regions is in good agreement with an estimate based on the force balance arguments. In the second part, we performed numerical simulations to study the effects of viscosity ratio, surface tension, contact angle, and drop size on the deformation, breakup, and critical pressure of permeation of oil droplets on membrane pores of circular cross-section. In our numerical setup, the oil droplet was exposed to a linear shear flow induced by the moving upper wall. The critical pressure of permeation was computed using a novel procedure where the total critical pressure was found as the sum of pressure jumps across the oilwater interfaces of the droplet in the channel and in the pore. First, the pressure jump across the static interface inside the pore was found using the Young-Laplace equation. Then, the pressure jump across the dynamic interface above the membrane surface was computed numerically and added to the pressure jump in the pore. This method was proven to be accurate, robust, and computationally efficient. To determine the dimensions of the 66 computational domain, we also studied the effect of confinement on the droplet deformation and breakup and concluded that in order to minimize finite size effects and computational costs, the distance between the membrane surface and the upper wall has to be at least twice the droplet diameter. In particular, it was observed that highly confined droplets become significantly deformed in a shear flow and break up more easily. In the absence of crossflow, we found that the analytical prediction for the critical permeation pressure derived by Nazzal and Wiesner [65] agrees well with the results of numerical simulations for different oil-to-water viscosity ratios, surface tension, contact angles, and droplet sizes. In general, with increasing crossflow shear rate, the critical permeation pressure increases with respect to its zero-shear-rate value and the droplet undergoes elongation in the flow direction followed by breakup into two segments. The results of numerical simulations indicate that at a fixed shear rate, the critical permeation pressure increases as a function of the viscosity ratio, which implies that more viscous droplets penetrate into the pore at higher transmembrane pressures. In agreement with a scaling relation for the critical capillary number, we also found that droplets of higher viscosity tend to break at lower shear rates. Furthermore, with increasing surface tension coefficient, the maximum increase in the critical permeation pressure due to crossflow becomes larger and the droplet breakup occurs at higher shear rates. Interestingly, the percent increase in critical permeation pressure as a function of the capillary number was found to be independent of the surface tension coefficient. Next, we showed that the breakup capillary number and the increase in critical pressure of permeation are nearly independent of the contact angle. Last, it was demonstrated that smaller droplets penetrate into the pore at lower pressures and break up at higher shear rates because larger shear stresses are needed to deform the interface above the membrane surface. 67 While most microfiltration membranes used in medium- to large-scale separation applications have pores of complex morphologies and a distribution of nominal sizes, results obtained for the simple case of a pore of circular cross-section can be useful for identifying general trends. As the model describes the interaction of oil droplets with unblocked pores, which corresponds to the initial stage of filtration or the stage that immediately follows membrane cleaning, the results can be helpful in understanding how fouling starts and in devising means to delay the onset of fouling. With the development of new methods of manufacturing micro-engineered membranes [125] and the rapid growth in the diversity and scale of applications of microfluidic devices, conclusions obtained in this work can be of direct practical value for guiding membrane design and optimizing process variables. 68 BIBLIOGRAPHY 69 BIBLIOGRAPHY [1] T.M. Squires, S.R. Quake, Rev. Mod. Phys. 77 (2005) 977. [2] J.C.T. Eijkel, A. van den Berg, Microfluid. Nanofluid. 1 (2005) 249. [3] M.J. Blunt, Curr. Opinion Colloid Interface Sci. 6 (2001) 197. [4] F.A. Coutelieris, J.M.P.Q. Delgado, Transport Processes in Porous Media, Springer, 2012. [5] L. Shui, J.C.T. Eijkel, A. van den Berg, Adv. Colloid Interface Sci. 133 (2007) 35. [6] S.-Y. Teh, R. Lin, L.-H. Hung, A.P. Lee, Lab Chip 8 (2008) 198. [7] C.D. McAuliffe, J. Petrol. Technol. 25 (1973) 721. [8] S.L. Kokal, B.B. Maini, R. Woo, Adv. Chem. 231 (1992) 220. [9] T. Sheorey, K. Muralidhar, P.P. Mukherjee, Int. J. Therm. Sci. 40 (2001) 981. [10] J. Mueller, Y. Cen, R.H. Davis, J. Membr. Sci. 129 (1997) 221. [11] M. Cheryan, N. Rajagopalan, J. Membr. Sci. 151 (1998) 13. [12] I.-S. Chang, C.-M. Chung, S.-H. Han, Desalination 133 (2001) 225. [13] G. Orive, R.M. Hernandez, A.R. Gascon, R. Calafiore, T.M.S. Chang, P. De Vos, G. Hortelano, D. Hunkeler, I. Lacik, A.M.J. Shapiro, J.L. Pedraz, Nature Medicine 9 (2003) 104. [14] J. Atencia, D.J. Beebe, Nature 437 (2005) 648. [15] J. Zhang, R.J. Coulston, S.T. Jones, J. Geng, O.A. Scherman, C. Abell, Science 335 (2012) 690. 70 [16] J. Ouellette, Indust. Phys. 9 (2003) 14. [17] H.A. Stone, A.D. Stroock, A. Ajdari, Annu. Rev. Fluid Mech. 36 (2004) 381. [18] D.J. Beebe, G.A. Mensing, G.M. Walker, Annu. Rev. Biomed. Eng. 4 (2002) 261. [19] D. Sparks, R. Smith, M. Straayer, J. Cripe, R. Schneider, A. Chimbayo, S. Anasari, N. Najafi, Lab Chip 3 (2003) 19. [20] H. Ohya, J.J. Kim, A. Chinen, M. Aihara, S.I. Semenova, Y. Negishi, O. Mori, M. Yasuda, J. Membr. Sci. 145 (1998) 1. [21] M.-F. Ficheux, L. Bonakdar, F. Leal-Calderon, J. Bibette, Langmuir 14 (1998) 2702. [22] L. L. Schramm, Adv. Chem. Series. 231 (1992). [23] L. C. Collins-Gold, R. T. Lyons, L. C. Bartholow, Adv. Drug Deliver. Rev. 5 (1990) 189. [24] B.T. Kelly, J.-C. Baret, V. Taly, A.D. Griffiths, Chem. Commun. 18 (2007) 1778. [25] D. Langevin, S. Poteau, I. Henaut, J.F. Argillier, Oil Gas Sci. Tech. 59 (2004) 511. [26] L.E. Sanchez, J.L. Zakin, Ind. Eng. Chem. Res. 33 (1994) 3256. [27] V. Schr¨oder, O. Behrend, H. Schubert, J. Coll. Int. Sci. 202 (1998) 334. [28] D. J. McClements, Crit. Rev. Food Sci. 47 (2007) 611. [29] J. Sjoblom, Surf. Sci. Series. 61 (1996). [30] R. Katoh, Y. Asano, A. Furuya, K. Sotoyama, M. Tomita, J. Membr. Sci., 113 (1996) 131. [31] M. Gryta, K. Karakulski, Desalination 121 (1999) 23. 71 [32] J. Bibette, F.L. Calderon, P. Poulin, Rep. Prog. Phys. 62 (1999) 969. [33] S. M. Joscelyne, G. Tr¨ag˚ ardh, J. Membr. Sci., 169 (2000) 107. [34] J. M. Montgomery, Wiley (1985). [35] B. Abisma¨ıl, J. P. Canselier, A. M. Wilhelm, H. Delmas, C. Gourdon, Ultrason. Sonochem. 6 (1999) 75. [36] H. Karbstein, H. Schubert, Chem. Eng. Process. 34 (1995) 205. [37] A. J. Gijsbertsen-Abrahamse, A. van der Padt, R. M. Boom, J. Membr. Sci. 230 (2004) 149. [38] P. Walstra, Chem. Eng. Sci. 48 (1993) 333. [39] G.T. Vladisavljevi`c, S. Tesch, H. Schubert, Chem. Eng. Process. 41 (2002) 231. [40] C.F. Christopher, S.L. Anna, J. Phys. D: Appl. Phys. 40 (2007) R319. [41] Y. Pan, W. Wang, T. Wang, P. Yao, Sep. Purif. Technol. 57 (2007) 388. [42] M. Gryta, K. Karakulski, A.W. Morawski, Water Res. 35 (2001) 3665. [43] L. Song, J. Membr. Sci. 139 (1998) 183. [44] J. Hermia, Trans. Inst. Chem. Eng. 60 (1982) 183. [45] P.M. Heertjes, Chem. Eng. Sci. 6 (1957) 190. [46] W.R. Bowen, J.I. Calvo, A. Hernandez, J. Membr. Sci. 101 (1995) 153. [47] V.T. Kuberkar, R.H. Davis, J. Membr. Sci. 168 (2000) 243. [48] J. Dufreche, M. Prat, P. Schmitz, Deslination 145 (2002) 129. [49] V.V. Tarabara, I. Koyuncu, M.R. Wiesner, J. Membr. Sci. 241 (2004) 65. 72 [50] R. W. Field, D. Wu, J. A. Howell, B. B. Gupta, J. Membr. Sci. 100 (1995) 259. [51] S.-H. Park, T. Yamaguchi, S. Nakao, Chem. Eng. Sci. 56 (2001) 3539. [52] A.J. Bromley, R.G. Holdich, I.W. Cumming, J. Membr. Sci. 196 (2002) 27. [53] A.B. Koltuniewicz, R.W. Field, T.C. Arnot, J. Membr. Sci. 102 (1995) 193. [54] R.G. Holdich, I.W. Cumming, I.D. Smith, J. Membr. Sci. 143 (1998) 263. [55] N. Bremond, J. Bibette, Soft Matter 8 (2012) 10549. [56] S. van der Graaf, T. Nisisako, C.G.P.H. Schroen, R.G.M. van der Sman, R.M. Boom, Langmuir 22 (2006) 4144. [57] J.H. Xu, G.S. Luo, G.G. Chen, J.D. Wang, J. Membr. Sci. 266 (2005) 121. [58] J. Husny, J.J. Cooper-White, J. Non-Newton. Fluid Mech. 137 (2006) 121. [59] A. Gupta, S.M. Sohel Murshed, R. Kumar, Appl. Phys. Lett. 94 (2009) 164107. [60] G. F. Christopher, N. N. Noharuddin, J. A. Taylor, S. L. Anna, Phys. Rev. E 78 (2008) 036317. [61] P. V. Puyvelde, A. Vananroye, R. Cardinaels, P. Moldenaers, Polymer 49 (2008) 5363. [62] H.A. Stone, Annu. Rev. Fluid Mech. 26 (1994) 65. [63] S. Lee, Y. Aurelle, H. Roques, J. Membr. Sci. 19 (1984) 23. [64] A. Hong, A. G. Fane, R. Burford, J. Membr. Sci. 222 (2003) 19. [65] F.F. Nazzal, M.R. Wiesner, Water Environ. Res. 68 (1996) 1187. [66] I.W. Cumming, R.G. Holdich, I.D. Smith, J. Membr. Sci. 169 (2000) 147. [67] G. I. Taylor, Proc. R. Soc. Lond. A 138 (1932) 41. 73 [68] G. I. Taylor, Proc. R. Soc. Lond. A 146 (1934) 501. [69] R. G. Cox, J. Fluid Mech. 37 (1969) 601. [70] E. J. Hinch, A. Acrivos, J. Fluid Mech. 98 (1980) 305. [71] F. D. Rumscheidt, S. G. Mason, J. Coll. Sci. 16 (1961) 238. [72] H.P. Grace, Chem. Eng. Commun. 14 (1982) 225. [73] B. J. Bentley, L. G. Leal, J. Fluid Mech. 167 (1986) 241. [74] J. M. Rallison, J. Fluid Mech. 109 (1981) 465. [75] C. Pozrikidis, Cambridge: Cambridge Univ. Press (1992). [76] J. U. Brackbill, D. B. Kothe, C. Zemach, J. Comp. Phys. 100 (1992) 335. [77] D. Gueyffier, J. Li, A. Nadim, R. Scardovelli, S. Zaleski, J. Comput. Phys. 152 (1999) 423. [78] J. Li, Y.Y. Renardy, M. Renardy, Phys. Fluids 12 (2000) 269. [79] V. Cristini, Y. Tan, Lab Chip 4 (2004) 257. [80] P. Dimitrakopoulos, J. Fluid Mech. 580 (2007) 451. [81] P. G. Saffman, J. Fluid Mech. 22 (1965) 385. [82] M. E. O’Neill, Chem. Eng. Sci. 23 (1968) 1293. [83] T. C. Price, Q. J. Mech. Appl. Math. 38 (1985) 93. [84] C. Pozrikidis, J. Eng. Math. 31 (1997) 29. [85] K. Sugiyama, M. Sbragaglia, J. Eng. Math. 62 (2008) 35. [86] A. Vananroye, P. Van Puyvelde, P. Moldenaers, Langmuir 22 (2006) 3972. 74 [87] H. Xi, C. Duncan, Phys. Rev. E 59 (1999) 3022. [88] A. Fakhari, M.H. Rahimian, Int. J. Numer. Meth. Fluids 64 (2010) 827. [89] X. Li, C. Pozrikidis, J. Fluid Mech. 341 (1997) 165. [90] P.J.A. Janssen, P.D. Anderson, Phys. Fluids 19 (2007) 043602. [91] L.J. Dietsche, A.C. Neubauer, Chem. Eng. Sci. 64 (2009) 4543. [92] V.R. Gopala, B.G.M. van Wachem, Chem. Eng. J. 141 (2008) 204. [93] J.E. Pilliod., E.G. Puckett, J. Comput. Phys. 199 (2004) 465. [94] Fluent, Inc., 2003. FLUENT 6.1 Users Guide. [95] C.W. Hirt, B.D. Nichols, J. Comput. Phys. 39 (1981) 201. [96] A.A. Saha, S.K. Mitra, J. Colloid Interface Sci. 339 (2009) 461. [97] N. Hsu, N. Ashgriz, J. Colloid Interface Sci. 270 (2004) 146. [98] D. Gerlach, G. Tomar, G. Biswas, F. Durst, Int. J. Heat Mass Transfer 49 (2006) 740. [99] W. J. Rider and D. B. Kothe, J. Comp. Phys. 141 (1998) 112. [100] M. Renardy, Y. Renardy, J. Li, J. Comput. Phys. 171 (2001) 243. [101] T. Darvishzadeh, N. V. Priezjev, J. Membr. Sci., 423-424 (2012) 468. [102] A. I. Hill, C. Pozrikidis, J. Colloid Interface Sci. 356 (2011) 763. [103] J. Eells, Math. Intell. 9 (1987) 53. [104] R. Finn, Capillary surface interfaces. Not. Am. Math. Soc. 46 (1999) 770. [105] P. Concus, R. Finn, Microgravity Sci. Technol. 3 (1990) 87. 75 [106] P. Concus, R. Finn, Acta Math. 132 (1974) 177. [107] V. Brady, P. Concus, R. Finn, Microgravity Sci. Technol. 15 (2004) 31. [108] A. Ullah, R. G. Holdich, M. Naeem, V. M. Starov, J. Membr. Sci. 401-402 (2012) 118. [109] N. Ichikawa, K. Hosokawa, R. Maeda, J. Colloid Interface Sci. 280 (2004) 155. [110] H. Cho, H.Y. Kim, J.Y. Kang, T.S. Kim, J. Colloid Interface Sci. 306 (2007) 379. [111] H. Wong, S. Morris, C.J. Radke, J. Colloid Interface Sci. 148 (1992) 317. [112] V.S. Ajaev, G.M. Homsy, Annu. Rev. Fluid Mech. 38 (2006) 277. [113] J. Feng, J.P. Rothstein, J. Colloid Interface Sci. 354 (2011) 386. [114] S.H. Collicott, M.M. Weislogel, AIAA 42 (2004) 289. [115] B.C. Berndt, Ramanujan’s Notebooks, Springer Verlag, 1985. [116] N. Albright, Report LBL-6137 (1977). [117] J.H. Xu, S.W. Li, J. Tan, G.S. Luo, Microfluid. Nanofluid. 5 (2008) 711. [118] E. van der Zwan, R. van der Sman, K. Schroen, R. Boom, J. Colloid Interface Sci. 335 (2009) 112. [119] Z. Zapryanov, S. Tabakova, Springer, Vol. 50 (1998). [120] Y. Renardy, Rheol. Acta. 46 (2007) 521. [121] J. M. Rallison, Annu. Rev. Fluid Mech. 16 (1984) 45. [122] H.A. Stone, PhD Thesis, Caltech (1988). ardh, C. Tra¨ag˚ ardh, J. Membr. Sci. 337 (2009) 232. [123] A. Timgren, G. Tra¨ag˚ 76 [124] K. Suzuki, I. Shuto, Y. Hagura, Food Sci. Tech. Int., 2(1) (1996) 43. [125] C.J.M. van Rijn, Micro-Engineered Membranes. In: Encyclopedia of Membrane Science and Technology, Eds: E.M.V. Hoek, V.V. Tarabara, John Wiley and Sons, (2013). [126] T. Darvishzadeh, V. V. Tarabara, N. V. Priezjev, J. Membr. Sci., 447 (2013) 442. 77