Numerical Simulations of Bubble Formation and Rise in Microchannels

Dec 31, 2008 - In this work, the formation of bubbles and their rise in circular ... The dynamics of formation and rise of Taylor bubbles in glass ...
0 downloads 0 Views 3MB Size
Download PDF
Ind. Eng. Chem. Res. 2009, 48, 8109–8120

8109

Numerical Simulations of Bubble Formation and Rise in Microchannels Deepak Goel and Vivek V. Buwa* Department of Chemical Engineering, Indian Institute of TechnologysDelhi, New Delhi 110 016, India

Gas-liquid flow in microchannels is of fundamental importance to many engineering applications involving microreactors, monolith reactors, microheat exchangers, and several other microfluidic devices. Slug flow, characterized by motion of long bubbles, also referred to as Taylor bubbles, is the most important of the different two-phase flow regimes observed in microchannels. In this work, the formation of bubbles and their rise in circular capillaries in the Taylor flow regime is investigated by using the volume-of-fluid method. The dynamics of formation and rise of Taylor bubbles in glass capillaries of 1, 0.5, 0.75, and 0.3 mm diameter for air-water and air-octane systems was simulated. The effects of superficial gas and liquid velocities, channel geometry (nozzle wall thickness, nozzle diameter, capillary diameter), wall adhesion (contact angle), and fluid properties (surface tension, viscosity) on the dynamics of bubble formation were investigated. The predicted bubble shapes and bubble formation periods were validated using the experimental data reported by Salman et al. (Salman, W.; Gavriilidis, A.; Angeli, P. Chem. Eng. Sci. 2006, 61, 6653-6666) for a wide range of experimental parameters. Such experimentally validated computational flow models will be useful to simulate the mass transfer and reactions in microcapillaries/channels. 1. Introduction Gas-liquid two-phase flow in microchannels is important to many process intensification applications through development and use of microreactors, monolith reactors, micromixers and heat exchangers, and several other microfluidic devices. Monolith catalysts, single blocks produced by extrusion of ceramic materials comprised of identical and parallel channels of different shapes (circular, rectangular, triangular), have been one of the important developments in catalysis. Monoliths offer several advantages over conventional reactors, particularly high specific surface area for mass/heat transfer and reactions, relatively lower pressure drops, lower axial dispersion and back mixing, and high catalyst efficiency because of the thin active layer of the catalyst. Developed in the 1970s for the automobile industry for the treatment of exhaust gases (NOx, CO from vehicles), monoliths have been tested for several demanding fast reactions at high temperatures such as steam reforming, partial oxidations, and oxidative dehydrogenations.2 For multiphase solid-catalyzed gas-liquid reactions, Nijhuis et al.3 have shown that monoliths outperform conventional reactors. Monoliths are already into commercial production, e.g., for peroxide synthesis and also for nitro-aromatic hydrogenations.4 With the progress in microfabrication technologies, the past decade witnessed significant research activity in the development of microreactors and the feasibility of their applications for small-scale productions, e.g., multichannel microreactors, micromixers, and micro heat exchangers. Microreactors provide high surface to volume ratios and therefore offer significant advantages of high heat and mass transport rates. Moreover, it is possible to carry out reactions under severe operating conditions and yet safely. Small reaction volumes can allow high temperature and concentration gradients, precise process control, and efficient heat management. As a result, a number of reaction systems have been tested using microreactor systems. These include several gas phase reactions, e.g., partial oxidations,hydrogenations,dehydrogenations,andreformingprocesses,5,6 where excellent temperature control and prevention of hot spots * To whom correspondence should be addressed. Tel.: +91 11 2659 1027. Fax: +91 11 2658 1120. E-mail: [email protected].

was achieved. As a consequence, increased selectivity toward desired products was obtained.6 Besides the single-phase reactions, various microreactor configurations have been evolved for gas-liquid, liquid-liquid, and gas-liquid-solid phase reactions, e.g., hydrogenations, chlorinations, and sulfonations.7 These include a falling film microreactor, multichannel and mesh microreactor, microbubble column, and micropacked bed reactor.7,8 While significant research is in progress to evolve different monolith structures and various configurations of microreactors and to test their feasibility for various multiphase reaction systems, the microscopic fluid dynamics which controls the interphase mass and heat transport and therefore the conversion and selectivity of an underlying reaction system is not yet well understood. With the development of efficient numerical methods and high-performance computing machines to solve basic equations of fluid flow, it is now possible to carry out detailed simulations of fluid flow and transport processes. While several numerical methods are being developed to simulate the multiphase flows with moving interfaces (see section 2.1), the capabilities of the numerical methods to predict the effects of gas and liquid superficial velocities, channel diameter, gas and liquid inlet configurations, surface characteristics of the channel wall, and physical properties of the liquid phase on the bubble formation dynamics and rise behavior in microchannels are not yet well investigated and experimentally verified. The present work is therefore focused on numerical simulations of gas-liquid flows in microchannels and their experimental validation. A brief review of the present state of art of experimental and numerical investigations of gas-liquid flows in microchannels is provided in section 2. The formulation of the volume-of-fluid (VOF) method used for the numerical simulations in the present work is discussed in section 3. The results are discussed in section 4, which is followed by the conclusions in section 5. 2. Previous Work Depending on the operating and design parameterssthe gas and liquid superficial velocities, the physical properties of fluid phases, the channel shape and size, the channel wall characteristics, and the gas and liquid inlet configurationsdifferent

10.1021/ie800806f CCC: $40.75  2009 American Chemical Society Published on Web 12/31/2008

8110

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

flow regimes, viz., bubbly, churn, Taylor, annular, and film flow, can be observed. The Taylor flow regime, characterized by long gas bubbles, sometimes also referred to as gas slugs, is one of the most important flow regimes observed under specific conditions. It has been shown that Taylor bubbles lead to enhancement in mass transfer because of good circulation within the liquid slugs and a thin characteristic liquid film between the bubble and walls.8 Over the past decade, a lot of work has been done on the experimental characterization of gas-liquid flows in microchannels by using different experimental techniques to understand the effects of the above-mentioned operating and design parameters on the flow regimes, key hydrodynamic parameters, e.g., bubble size/shape, their rise velocity, film thickness, slug lengths, pressure drop, residence time distribution, and gas-liquid mass transfer. A brief review of the above-mentioned experimental investigations and also of numerical investigations, relevant to present study, is provided below. 2.1. Experimental Investigations. Unlike the conventional columns/channels, invasive measurement techniques such as intrusive probes have limited applicability to microchannels because of their small size. A number of studies in microchannels have been carried out experimentally using different noninvasive techniques such as the following: γ-ray computed tomography (CT) has been used for studying gas/liquid holdup distribution in monoliths;10 magnetic resonance imaging (MRI) has been used for quantitative visualization of gas-liquid flow;11,12 NMR microimaging was found suitable for studying flow characteristics in complicated shaped channels;13 microparticle image velocimetry (µ-PIV) provided velocity field distribution inside liquid slugs;14-16 flow visualization used high-speed imaging to determine flow regime maps, slug lengths, liquid film thickness, average void fraction and bubble formation mechanisms,17-20 and residence time distributions (RTDs).21 Several investigators studied flow regimes and reported the flow regime maps for a wide range of gas and liquid velocities (0.01 < UG, UL < 10 m/s) and for different channel geometries (circular, rectangular, triangular cross-section channels, 50 µm < DH < 3 mm).16,19,22-28 In most of these investigations, flow regime maps were studied as a function of superficial gas and liquid velocities or flow rates for a particular gas-liquid system and channel geometry. Therefore, the generalized criteria of the regime transition characteristics with the change in UG, UL, or DH cannot be established from the reported literature. Waelchli et al.16 presented generalized flow regime maps based on Reynolds number (Re ) DHUF/µ) and Weber number (We ) DH(FL - FG)U2/σ), taking into account the effects of superficial velocities, the density of the fluid, the surface tension, the viscosity, and the hydraulic diameter of the channel. However, the effects of channel cross section and surface characteristics (contact angle and roughness) still remain unaccounted for. Coleman and Garimella29 experimentally investigated gas-liquid flow in channels with circular and rectangular cross sections (1.3 mm < DH < 5.5 mm) and reported a strong dependency of flow regimes on the shape and diameter of the channel. Therefore, the utility of the flow regime maps for a priori identification of the flow regimes for different microchannels (other than the channels for which the map was experimentally obtained) is highly dependent on the similarity of both the channel cross-sectional shape and dimensions rather than only on a similarity in hydraulic diameter. The Taylor flow regime is characterized by the capsule-shaped long bubbles separated by liquid slugs and a thin liquid film

between the bubbles and the channel wall. The gas and liquid slug lengths strongly depend on the operating and design parameters mentioned above. Kreutzer et al.30 suggested the following correlation for estimating the slug length using the pressure drop measurements (∆P/L): fRe )

[

16 0.17d Re 1+ Re Lslug Ca

( )

1/3

]

(1)

where d is the diameter of the channel, Lslug is the liquid slug length, and fRe is the friction factor at a particular Re. Here it must be noted that fRe in itself depends on the value of ratio Lslug/d. They measured slug lengths for 0.02 m/s < UL < 0.2 m/s and 0.02 m/s < UG < 0.3 m/s in a channel of 1.56 mm diameter for different types of gas-liquid inlet distributors at the inlet, and they observed that slug lengths also vary significantly with the type of distributor and hence cannot be predicted accurately by using the correlations that do not account for the influence of inlet gas and liquid distribution on the bubble formation process. The pressure drop inside the channel and the thickness of the thin liquid film between the bubble and channel wall are two other important flow variables required for the design of microchannels. Woehl and Cerro31 proposed a theoretical model for the computation of pressure drop in capillaries for the Taylor flow regime based on hydrostatics pressure drop, viscous pressure drop, and capillary pressure drop. Mishima and Hibiki24 investigated a two-phase frictional pressure drop as a function of the inner diameter of the channel. Pehlivan et al.32 measured a two-phase frictional pressure drop in channels of 3, 1, and 0.8 mm diameters and observed that, as the channel diameter decreases, existing pressure drop models for conventional reactors become inaccurate for predicting the pressure drop inside the microchannels. A very thin liquid film (5-50 µm) between the bubble and the channel wall provides a region of high mass transfer rates from the gas to the wall because of the large interfacial area, short diffusion path, and high liquid velocities inside the region. The liquid film thickness was studied experimentally by Taylor,33 Marchessault and Mason,34 Brethreton,35 Goldsmith and Mason,36 and Schwart et al.37 Waelchili et al.15 analyzed liquid recirculation and presented a detailed study of the liquid velocity distribution in intermittent flows in microchannels. van Stejin et al.14 studied the bubble formation mechanism and provided a detailed velocity flow field using µ-PIV. Salman et al.1 studied different mechanisms of the formation of the Taylor bubbles and reported the bubble formation periods for different inlet nozzles and for a wide range of gas and liquid superficial velocities. Several studies have been done on understanding various parameters influencing the heat and mass transfer phenomenon in capillaries.38-41 It has been shown that the heat and mass transfer characteristics of the two-phase flow in capillaries largely depends on the gas-liquid slug lengths, liquid film thickness, bubble shape, and bubble velocity among other parameters. For instance, Oliver and Wright42 and Horvath et al.43 demonstrated a radial increase in mass and heat transfer for the slug flow inside capillaries. Bercˇic and Pintar44 studied the role of gas-liquid slug lengths in mass transfer and reported that the major contribution of the mass transfer is achieved through the surface exposed to liquid rather than through gas bubbles. 2.2. Numerical Investigations. Several methods have been developed to simulate gas-liquid flows with moving interfaces, mainly, the front tracking (FT) method,45,46 the level set (LS) method,47 the VOF method,48-50 and the lattice Boltzmann (LB)

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009 51,52

method. It should be noted that the above-mentioned methods have been largely applied to simulate the free rise of single/multiple bubbles and predictions have been experimentally validated. However, the predictive capabilities of these methods to simulate the gas-liquid flows in microchannels have not been fully explored. A brief review of the applications of the above-mentioned numerical methods and the present state of the art of their predictive capabilities to simulate the gas-liquid flows in microchannels is provided below. Over the past few years, only a few numerical investigations have been reported on numerical computations of gas-liquid flow in microchannels.9,53-60 Primarily, two approaches were adapted for the simulation of gas-liquid flow in microchannels. In the first approach, the initial shape of a bubble was defined, and the frame of reference was fixed with respect to the bubble; i.e., the bubble was stationary and the wall was moving with the velocity of the bubble.53-58 The focus was on the simulation of shape evolution of a single bubble with presumed bubble length. van Baten and Krishna56,57 performed axisymmetric VOF simulations to investigate the effects of bubble velocity, liquid film thickness, and slug length on the mass transfer from a bubble to the liquid phase and from the liquid phase to the channel wall for the Taylor flow regime. Taha and Cui,53-55 using two-dimensional (2-D) VOF simulations, investigated the effects of channel geometry on bubble shape, bubble rise velocity, and wall shear distribution in circular and square capillaries. Yang et al.,58 using the 2-D LB method, investigated the flow regimes and the effects of body force, surface tension, and bubble size on the rise behavior of Taylor bubbles in microchannels. Irandoust and Anderson8 simulated Taylor flow in capillaries using the finite difference method and demonstrated the enhancement of mass transfer because of the recirculation within the characteristic thin liquid film and the liquid slug. In the second approach, the dynamics of bubble formation and rise of bubbles in microchannels was simulated.59,60 Qian and Lawal59 performed 2-D VOF simulations of bubble formation and their rise in square channels with cross-sectional widths of 0.25, 0.5, 0.75, 1, 2, and 3 mm and compared the predicted gas and liquid slug lengths using the experimental data of Kreutzer et al.30 The predicted gas and liquid slug lengths were found to be in qualitative agreement with the measurements. However, the effects of the gas-liquid-wall contact angle and the exact geometry of the gas and liquid inlet as used by Kreutzer et al.30 were not accounted for. Qian and Lawal59 further showed that the predicted slug lengths were highly sensitive to the gas and liquid inlet configuration, besides the channel diameter and shape. They observed that the simulated Laplace pressure for a three-dimensional (3-D) domain was twice that obtained for a 2-D domain, because of the curvature effects. Qian and Lawal59 also reported that the predictions of the 2-D and 3-D computations were comparable and that the 2-D simulations were adequate. However, as shown later in this work, at coarse grid resolutions, the predictions of both 2-D and 3-D computations are comparable and agree well with the experiments. However, such an agreement observed at coarse grid resolution worsens with further grid refinements, and predictions of 2-D simulations are not grid independent. Therefore, 3-D effects must be considered to predict bubble formation mechanism and rise behavior accurately. Yu et al.60 performed 2-D LB simulations to investigate bubble formation dynamics, bubble shape, and bubble size for different geometries and observed that gas-liquid inlet configuration has a significant influence on the bubble formation mechanism and flow char-

8111

acteristics. The simulation results were in reasonable qualitative agreement with the experiments, and the observed differences in the predictions and experiments were attributed to the inaccurate representation of the wetting conditions and the 3-D effects. The role of wall adhesion in flow characteristics such as slug lengths and film thickness has been broadly accepted, but was not considered in the previous studies. To the best of the authors’ knowledge, in most of the previous numerical investigations the effects of wall adhesion were ignored. Since the gas-liquid flow in microchannels strongly depends on the inlet configuration, the exact inlet geometry of the gas and liquid inlet needs to be considered and wall characteristics need to be precisely accounted for the purpose of rigorous experimental validation. The objective of the present work is to investigate the bubble formation dynamics and the rise behavior of bubbles in microchannels. In the present work, we have numerically investigated the formation of Taylor bubbles in glass capillaries of 1, 0.75, 0.5, and 0.3 mm diameter using the VOF method. The numerical simulations were carried out using 2-D axisymmetric and 3-D solution domains, and the predictions were compared with the bubble shapes and formation periods reported by Salman et al.1 The effects of wall surface adhesion (contact angle), inlet configuration (nozzle diameter, nozzle wall thickness and contact angle), gas superficial velocity (0.008-0.03 m/s), liquid superficial velocity (0.0025-0.025 m/s), and liquid phase properties (viscosity 0.0001-0.01 Pa · s and surface tension 0.007-0.1 N/m) are investigated in detail. 3. Computational Model In this study, the VOF method based on the piecewise linear interface representation (PLIC-VOF) was used. The VOF method, originally proposed by Hirt and Nichols,48 is an efficient numerical simulation technique for complex free surface problems. Gas and liquid phases were modeled as incompressible fluids with constant values of viscosity and surface tension, and the flow was assumed to be Newtonian and laminar. 3.1. The VOF Method. In the VOF method, a particular fluid phase is defined by the volume fraction (Rq) of the qth phase inside the cell. If Rq is 0, the cell is empty of the qth phase, and if Rq is 1, the cell is full of the qth phase. The cells with 0 < Rq < 1 contain the interface between the qth phase and the other (pth) phase. Depending upon Rq, a single set of mass and momentum equations (eqs 2 and 3, respectively) for an incompressible Newtonian fluid and laminar flow was solved in the entire computational domain. ∂ (F) + ∇·(Fυb) ) 0 ∂t

(2)

∂ (Fυb) + ∇·(Fυb υb) ) -∇p + µ(∇2υb ) + Fg b+b F (3) ∂t where b υ is the velocity vector, p is the pressure, and b F is the surface tension force per unit volume. When a computational cell contains the gas-liquid interface, the volume fraction weighted average properties of the gas and liquid phases were used as F ) RGFG + (1 - RG)FL

(4)

µ ) RGµG + (1 - RG)µL

(5)

The surface tension force was modeled using the continuum surface force (CSF) model proposed by Brackbill et al.49 With this model, the surface tension force results in a source term in the momentum equation. The interface curvature κ was computed as

8112

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

κ ) -(∇·nˆ) (6) where nˆ is the unit normal and calculated as nˆ ) n/|n| and n ) b) in eq 3 was ∇Rq. The surface tension force per unit volume (F calculated as b F ) σpq

Fκq∇Rq (1/2)(Fp + Fq)

(7)

where F is the volume-averaged density, σ is the surface tension coefficient, and κ is the interface curvature. The subscripts p and q refer to the particular fluid phase. The position of the gas-liquid interface was advanced by solving an advection equation for Rq as48 ∂Rq + υb ·∇Rq ) 0 (8) ∂t The geometric reconstruction scheme based on the piecewise linear interface construction (PLIC)50,61 was used. In the 2-D PLIC method, it is assumed that the interface has linear slope within each cell and can be represented by the line nˆ·x + φ ) 0 (9) in which nˆ is the unit normal vector, x is a point on the line, and φ is the line constant. The values of nˆ and φ were calculated for each cell containing the interface. For fluids exhibiting a nonzero contact angle, the presence of the wall also affects the surface normal and thus wall adhesion forces were included in the present work. The three-phase contact angle was used to adjust the surface normal in the cells near the wall as nˆ ) nˆw cos θw + ˆtw sin θw

(10) ˆ where nˆw and tw are the unit vectors normal and tangential to the wall, respectively. It was assumed that the contact line maintains a constant (static) contact angle (θw) independent of the velocity and the direction of the contact line movement. This boundary condition adjusts the curvature of the interface near the wall, which is in turn used in the surface tension force calculation. 3.2. Numerics, Solution Domain, and Boundary Conditions. The model equations (2), (3), and (8) given above were solved using the commercial flow solver Fluent 6.2. In the present work, circular microchannels of 0.3, 0.5, 0.75, and 1 mm diameters were used to investigate the formation of bubbles and their rise behavior. Simulations were carried out using both the 2-D axisymmetric and 3-D solution domains as shown in Figure 1. A no-slip boundary condition was used at the channel and nozzle wall(s), the velocity inlet boundary condition was used at the gas and liquid inlets and the inlet gas and liquid velocities were specified, the outlet of the channel was specified as outflow, which treats the flow as fully developed at the exit, and the static gas-liquid-wall three-phase contact angle was specified at the channel wall (θwall) and at the nozzle wall (θnozzle). Later, it was shown that, before the final detachment, the bubble base advances on and then recedes from the nozzle wall. A no-slip boundary condition was imposed for the fluid velocities at the nozzle and capillary wall. However, the moving contact line feels a slip. In order to avoid the stress singularities at the contact point, often a Navier slip condition is used.62,63 In most applications, this slip length is expected to be much smaller than the mesh size used in the numerical computations.62 In the present work, the available computational resources do not allow the mesh size to be refined in such a way that the true slip length is captured correctly. In the present work, therefore, no explicit slip length was imposed. However, it should be noted that the algorithm itself introduces an implicit

Figure 1. (a) Axisymmetric solution domain. (b) Three-dimensional solution domain and boundary conditions.

slip length on the order of the grid spacing (see Renardy et al.62 and Gerlach et al.64 for further details). The experimental device used by Salman et al.1 consisted of a 30 cm long capillary and a 4 cm long nozzle to ensure fully developed flow of the gas phase and the liquid phase. Axisymmetric simulations were done to study the effect of the inlet velocity profile (uniform and parabolic velocity profile) and the difference in predicted bubble formation periods was within 10 ms, which was less than 10% for the cases investigated in this work. Initial simulations were performed to study the effect of the discretization scheme (second-order UPWIND and QUICK) and pressure-velocity coupling algorithms (SIMPLE and PISO) on the predicted bubble formation period and bubble shapes. In both cases, the difference in the predicted bubble formation period was less than 10%. However, in the case of QUICK and PISO, the gas-liquid interface was sharper than second-order UPWIND and SIMPLE, respectively. However, the computational times for PISO were significantly higher than those with SIMPLE. Therefore, in all the simulations the QUICK scheme was used for the discretization of the momentum equation, the first-order implicit method was used for the discretization of the temporal derivatives, the second-order scheme was used for the pressure interpolation, the SIMPLE algorithm was used for the pressure-velocity coupling, and the geometric reconstruction scheme was used for the interface reconstruction. Depending on the grid resolution, and UG and UL, a time step varying between 1 × 10-6 and 1 × 10-7 s was used in the simulations. The adequacy of 2-D axisymmetric and 3-D solution domains and the effects of grid resolution were investigated, and the results are discussed separately in section 4.1. 4. Results and Discussion Salman et al.1 experimentally investigated the gas-liquid flow in a glass capillary with an internal diameter (dcap) of 1 mm for which the gas was introduced through nozzles of internal diameters (dnozzle) of 0.34 and 0.11 mm with outside diameters (o.d.’s) of 0.64 and 0.21 mm, respectively. The liquid was introduced through the annular space between the nozzle wall and capillary wall (see Figure 1). Salman et al.1 investigated

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

the effects of superficial gas velocity (0.002-0.04 m/s) and superficial liquid velocity (0.001-0.025 m/s) for air-water and air-octane systems (water: FL ) 1000 kg/m3, µL ) 8.5 × 10-4 Pa · s, σ ) 0.072 N/m; octane: FL ) 703 kg/m3, µL ) 0.001 Pa · s, σ ) 0.021 N/m; air: FG ) 1.19 kg/m3, µG ) 1.8 × 10-5 Pa · s) and reported bubble shapes and bubble formation periods. The contact angles for air-water-glass (θwall ) 10°) and air-water-steel (θnozzle ) 80°, steel nozzle was used in the experiments)wereexperimentallymeasured.Forair-octane-glass and air-octane-steel systems, octane was found to completely wet both materials and the contact angle was found to be very small (θwall ∼ θnozzle ∼ 1°). 4.1. Effect of Nozzle Wall Thickness and Solution Domain. The preliminary VOF simulations were carried out to investigate the importance of the nozzle geometry, the adequacy of the axisymmetric solution domain, and the effect of grid resolution. In the preliminary simulations, a standard test case of gas-liquid flow through a capillary (dcap ) 1 mm and lcap ) 5 mm), in which the gas (air) was introduced through a nozzle with dnozzle ) 0.34 mm (o.d. ) 0.64 mm) and with UG ) 0.012 m/s and the liquid (octane) was introduced with UL ) 0.0178 m/s, was considered and the bubble shapes and formation periods reported by Salman et al.1 were used to validate the predictions. In the simulations carried out without considering the thickness of the nozzle wall, i.e., gas entered the channel through a circular nozzle of zero wall thickness, it was observed that the bubble shapes were highly unstable and were detached immediately after 2-5 ms whereas the corresponding experimental bubble formation period was 76 ms. However, when the inlet geometry was modified by considering the thickness of the nozzle wall, the predicted bubble shapes were found to be stable and well developed with the realistic bubble formation period of 83 ms, as shown in Figure 2a. This indicates that the dynamics of bubble spreading and receding on the nozzle wall is important and should be considered in all further simulations. Though the predicted bubble formation period of 83 ms (obtained from axisymmetric simulation with 8500 cells corresponding to ∼34 cells per capillary diameter (cpd)) was in satisfactory agreement with the measured formation period of 76 ms (Salman et al.,1 as shown in Figure 2b), it is necessary to check the grid dependency of the predictions. Further simulations were carried out with 8500 (∼34 cpd), 34 000 (∼68 cpd), 77 500 (∼96 cpd), and 136 000 (∼135 cpd) cells, and the predicted formation periods of the first bubble were found to be 83, 113, 180, and 250 ms, respectively, against the measured formation period of 76 ms. One possible reason for the significant variation in the predicted bubble formation periods with the grid resolution could be that the curvature effects are not correctly simulated. In the process of bubble formation, the gas bubble was found first to spread on the nozzle wall, to grow vertically, and then to recede, before the neck formation and the final detachment occurred. The neck shapes observed just before the detachment for the above-mentioned grid resolutions (shown in Figure 3) indicate that, as the grid was refined, the neck length was also found to increase and the formation periods were also proportionally increased. Qian and Lawal59 compared the Laplace pressure predicted by both the 2-D and 3-D simulations and found that the Laplace pressure predicted by the 3-D simulations was twice that of the 2-D simulations, and attributed it to the curvature effects. Therefore, the axisymmetric simulations are not adequate to investigate the bubble formation mechanism in microchannels. Yu et al.,60 who investigated the gas-liquid flow in microchannels using 2-D lattice Boltzmann computations, also attributed the dis-

8113

Figure 2. Comparison of (a) predicted bubble shapes and formation period using axisymmetric solution domain with (b) the experiments of Salman et al.1 [reprinted with permission from ref 1, copyright Elsevier, 2006] for air-octane system (dcap ) 1 mm, lcap ) 5 mm, dnozzle ) 0.34 mm, UG ) 0.012 m/s, UL ) 0.0178 m/s, θnozzle ) 70°, θcap ) 40°, grid ) 8500 cells).

Figure 3. Inadequacy of axisymmetric simulations. Effect of grid resolution on bubble formation process for air-octane system (dcap ) 1 mm, lcap ) 5 mm, dnozzle ) 0.34 mm, UG ) 0.012 m/s, UL ) 0.0178 m/s, θnozzle ) 70°, θcap ) 40°): (a) 8500, (b) 34 000, (c) 77 500, and (d) 136 000 cells.

crepancy in the predicted and experimental bubble shapes to the inadequacy of 2-D simulations. Further, it should also be noted from the experiments of Salman et al.1 that the bubble formation is not exactly symmetric and a clear inclination of

8114

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

Figure 4. Predicted bubble shapes and formation period obtained using the 3-D solution domain [see for comparison the bubble shapes observed experimentally by Salman et al.1 shown in Figure 2b] for air-octane system (dcap ) 1 mm, dnozzle ) 0.34 mm, lcap ) 5 mm, UG ) 0.012 m/s, UL ) 0.0178 m/s, θnozzle ) 70°, θcap ) 40°, grid ) 33 558 cells).

Figure 6. Effect of superficial gas and liquid velocities on predicted bubble shapes for air-water system: (a) UL ) 0.0178 m/s and (i) UG ) 0.0120 m/s, (ii) UG ) 0.0154 m/s, and (iii) UG ) 0.0284 m/s; (b) UG ) 0.012 m/s and (i) UL ) 0.0075 m/s, (ii) UL ) 0.0178 m/s, and (iii) UL ) 0.023 m/s (dcap ) 1 mm, dnozzle ) 0.34 mm, lcap ) 5 mm, θnozzle ) 80°, θcap ) 10°). Figure 5. Adequacy of 3-D solution domain. Effect of grid resolution on the bubble formation period for air-octane system: (a) 33 558 cells (bubble formation period ) 76 ms), (b) 69 564 cells (bubble formation period ) 79 ms), and (c) 151 260 cells (bubble formation period ) 77 ms) (dcap ) 1 mm, dnozzle ) 0.34 mm, lcap ) 5 mm, UG ) 0.012 m/s, UL ) 0.0178 m/s, θnozzle ) 70°, θcap ) 40°).

the bubble toward the right-side channel wall is visible at 40 ms (see Figure 2b). In order to investigate further, 3-D simulations of bubble formation were carried out at different grid resolutions, corresponding to the total of 33 558 cells (∼18 cells per diameter and ∼102 cells along the length), 69 564 cells (∼29 cells per diameter and ∼102 cells along the length), and 151 260 cells (∼33 cells per diameter and ∼170 cells along the length). The predicted bubble shapes obtained from the simulation carried out with 33 558 grid cells are shown in Figure 4, which agreed very well with the experimental images of Salman et al.1 (see Figure 2b). In the 3-D simulations the asymmetric growth of the bubble can also be seen in Figure 4, which agreed well with the experiments (see Figure 2b). Importantly, the predicted formation periods for the simulations carried out with total cells of 33 558, 69 564, and 151 260 were 76, 79, and 77 ms, respectively (see Figure 5), and indicated that predictions were independent of the grid resolution (unlike the grid-dependent detachment behavior observed in the axisymmetric simulations). It should be also noted that the maximum frame grabbing speed

of 250 frames/s was used in the experiments and therefore a minimum measurement error of (4 ms in the formation periods is expected. In order to fix the length of the solution domain, simulations were done to study the effect of channel length on the bubble formation period. The bubble formation periods for the simulations carried out with lcap ) 5, 10, and 20 mm were found to be 76, 79, and 78 ms, respectively, and hence indicated no significant change in the bubble formation period with the channel length. Hence, to reduce computational time, all further simulations were carried out with lcap ) 5 mm (lcap/dcap ) 5). In addition, it was observed that bubble formation periods vary slightly for the successive bubble detachments, which is indicated by the error bar in the experimental data reported by Salman et al.1 (see Figures 7 and 8 below). The variation between the bubble formation periods widens as the liquid velocity is decreased. Even in simulations, the successive bubble formation periods were also found to vary within (4 ms. Qian and Lawal59 also reported a distribution in gas slug lengths and observed that slug length variation became wider as the superficial gas or liquid velocity was increased. 4.2. Effect of Gas and Liquid Superficial Velocities. The influence of gas and liquid superficial velocities on the bubble formation period was investigated for air-water and air-octane systems. The effect of UG and UL was first investigated for the air-water system (dcap ) 1 mm, dnozzle ) 0.34 mm, lcap ) 5 mm, θnozzle ) 80°, θcap ) 10°). The predicted bubble shapes

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

8115

Figure 7. Comparison of predicted bubble formation periods for air-water system for nozzle of dnozzle ) 0.34 mm with the experiments of Salman et al.1 (dcap ) 1 mm, lcap ) 5 mm, θnozzle ) 80°, θcap ) 10°).

Figure 9. Comparison of (a) predicted bubble shapes and formation period for air-water system for nozzle of dnozzle ) 0.11 mm with (b) experiments of Salman et al.1 [reprinted with permission from ref 1, copyright Elsevier, 2006] for air-water system (dcap ) 1 mm, lcap ) 5 mm, UG ) 0.012 m/s, UL ) 0.0178 m/s, θnozzle ) 80°, θcap ) 10°).

Figure 8. Comparison of predicted bubble formation periods for air-octane system for nozzle of dnozzle ) 0.34 mm with the experiments of Salman et al.1 (dcap ) 1 mm, lcap ) 5 mm, θnozzle ) 1°, θcap ) 1°).

obtained for different UG values (and at UL ) 0.0178 m/s) and for different UL values (and at UG ) 0.012 m/s) are shown in Figure 6a and 6b, respectively. A quantitative comparison of predicted and measured1 bubble formation periods is shown in Figure 7. For most cases (UL > 0.012 m/s), the average deviation of numerically predicted bubble formation periods was less than 10% of the experimental data. Figure 7 shows that the predicted bubble formation period decreases with increases in UG and UL. At lower values of UL (0.012 m/s) and low UG ( 0.012 m/s and UG < 0.02 m/s. However, for UL < 0.0075 m/s and UG > 0.02 m/s, bubble modification due to bubble pairing and coalescence at the nozzle was observed and deviation of numerical predictions from experimental data increased, though the experimental data contain large deviations (as large as 50%) in such cases.

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

8119

Notation d ) diameter, m DH ) hydraulic diameter of the channel, m b F ) continuum surface force, kg · m-2 · s-2 f ) friction factor g ) gravitational constant, m/s2 l ) length of the channel, m Lslug ) liquid slug length, m nˆ ) unit normal vector P ) pressure, Pa q ) particular phase under consideration t ) time, s ˆt ) unit tangent vector U ) superficial velocity, m/s We ) Weber number ()FdU2/σ) Re ) Reynolds number ()FdU/µ) Ca ) capillary number ()µU/σ) Greek Symbols R ) volume fraction of the phase in the computational cell κ ) surface curvature µ ) molecular viscosity, Pa · s θ ) contact angle, deg F ) density, kg/m3 σ ) surface tension, N/m b υ ) velocity vector, m/s φ ) line constant Subscripts and Superscripts G ) gas L ) liquid cap ) capillary w ) wall Figure 18. Effect of liquid phase physical properties on bubble formation period for air-water system: (a) surface tension and (b) viscosity (dcap ) 1 mm, dnozzle ) 0.34 mm, lcap ) 5 mm, UG ) 0.012 m/s, UL ) 0.0178 m/s, θnozzle ) 80°, θcap ) 10°).

(d) The wall adhesion was found to strongly influence the bubble formation mechanism and must be considered in the simulations for accurate prediction of bubble formation dynamics. The bubble formation period was found to decrease with the increase in the capillary wall contact angle and to increase with the increase in the nozzle wall contact angle. (e) The predicted bubble formation periods were found to strongly depend on the surface tension up to 0.072 N/m. With the increase in the surface tension, the bubble tends to become more spherical; the bubble diameter and hence the bubble formation period were found to increase. For the range of liquid viscosities considered in the present work (0.0001-0.01 Pa · s), the effect of viscosity on bubble formation was found to be negligible, which agreed well with the previous investigations reported in the literature. (f) Bubble formation is sensitive to gas and liquid inlet geometry and channel dimensions. The bubble formation period was found to decrease with the decrease in the nozzle diameter and the capillary diameter. The extensive validation of the computed bubble shapes and formation periods with the experimental results demonstrate that the Taylor flow regime in microchannels can be successfully predicted for a wide range of operating parameters. The experimentally validated computational model will be very useful to simulate the transport processes and reactions in microcapillaries/channels.

Literature Cited (1) Salman, W.; Gavriilidis, A.; Angeli, P. On the formation of Taylor bubbles in small tubes. Chem. Eng. Sci. 2006, 61, 6653–6666. (2) Kapteijn, K.; Nijhuis, T. A.; Heiszwolf, J. J.; Moulijn, J. A. New non-traditional multiphase catalytic reactors based on monolithic structures. Catal. Today 2001, 66, 133–144. (3) Nijhuis, T. A.; Kreutzer, M. T.; Romijn, A. C. J.; Kapteijn, F.; Moulijn, J. A. Monolithic catalysts as more efficient three-phase reactors. Catal. Today 2001, 66, 157–165. (4) Kreutzer, M. T.; Kapteijn, F.; Moulijn, J. A. Shouldn’t catalysts shape up? Structured reactors in general and gas-liquid monolith reactors in particular. Catal. Today 2006, 111, 111–118. (5) Lowe, H.; Ehrfeld, W. State-of-the-art in microreaction technology: concepts, manufacturing and applications. Electrochim Acta 1999, 44, 3679-3689. (6) Kolb, G.; Hessel, V. Micro-structured reactors for gas phase reactions. Chem. Eng. J. 2004, 98, 1–38. (7) Hessel, V.; Lo¨b, P.; Lo¨we, H. Industrial and real-life applications of micro-reactor process engineering for fine and functional chemistry. Stud. Surf. Sci. Catal. 2006, 159, 35–46. (8) Jensen, K. F. Microreaction engineeringsis small better? Chem. Eng. Sci. 2001, 56, 293–303. (9) Irandoust, S.; Andersson, B. Simulation of flow and mass transfer in Taylor flow through a capillary. Comput. Chem. Eng. 1989, 13, 519– 526. (10) Roy, S.; Al-Dahhan, M. Flow distribution characteristics of a gasliquid monolith reactor. Catal. Today 2005, 105, 396–400. (11) Gladden, L. F.; Lima, M. H. M.; Mantle, M. D.; Sederman, A. J.; Stitt, E. H. MRI visualization of two-phase flow in structured supports and trickle-bed reactors. Catal. Today 2003, 79-80, 203–210. (12) Sederman, A. J.; Heras, J. J.; Mantle, M. D.; Gladden, L. F. MRI strategies two phase flow in parallel channel ceramic monoliths. Catal. Today 2007, 128, 3–12. (13) Koptyug, I. V.; Ilyina, L. U.; Matveev, A. V.; Sagdeev, R. Z.; Parmon, V. N.; Altobelli, S. A. Liquid and gas flow and related phenomena

8120

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

in monolithic catalysts studied by H NMR micro imaging. Catal. Today 2001, 69, 385–392. (14) Tsoligas, A. N.; Simmons, M. J. H.; Wood, J. Two phase gasliquid reaction studies in a circular capillary. Chem. Eng. Sci. 2007, 62, 4365–4378. (15) van Steijn, V.; Kreutzer, M. T.; Kleijn, C. R. µ-PIV study of the formation of segmented flow in micro fluidic T-junctions. Chem. Eng. Sci. 2007, 62, 7505–7514. (16) Waelchli, S.; von Rohr, P. R. Two-phase flow characteristics in gas-liquid micro reactors. Int. J. Multiphase Flow 2006, 32, 791–806. (17) Serizawa, A.; Feng, Z.; Kawahara, Z. Gas-liquid two-phase flow in micro channels Part I: Two-phase flow patterns. Exp. Therm. Fluid Sci. 2002, 26, 703–714. (18) Kawahara, Z.; Chung, P. M.-Y.; Kawaji, M. Investigation of twophase flow pattern, void fraction and pressure drop in a micro channel. Int. J. Multiphase Flow 2002, 28, 1411–1435. (19) Hetsroni, G.; Mosyak, A.; Segal, Z.; Pogrebnyak, E. Two-phase flow patterns in parallel micro-channels. Int. J. Multiphase Flow 2003, 29, 341–360. (20) van Steijn, V.; Kreutzer, M. T.; Kleijn, C. R. Velocity fluctuations of segmented flow in microchannels. Chem. Eng. J. 2008, 135, 159–165. (21) Bakker, J. J. W.; Kreutzer, M. T.; de Lathouder, K.; Kapteijn, F.; Moulijn, J. A.; Wallin, S. A. Hydrodynamic properties of a novel ‘open wall’ monolith reactor. Catal. Today 2005, 105, 385–390. (22) Fukano, T.; Kariyasaki, A. Characteristics of gas-liquid two-phase flow in a capillary tube. Nucl. Eng. Des. 1993, 141, 59–68. (23) Xu, J. L.; Cheng, P.; Zhao, T. S. Gas-liquid two-phase flow regimes in rectangular channels with mini/micro gaps. Int. J. Multiphase Flow 1999, 25, 411–432. (24) Mishima, K.; Hibiki, T. Some characteristics of air-water two phase flow in small diameter vertical tubes. Int. J. Multiphase Flow 1996, 22 (4), 703–712. (25) Triplett, K. A.; Ghiaasiaan, S. M.; Abdel-Khalik, S. I.; Sadowski, D. L. Gas liquid two-phase flow in microchannels Part I: two-phase flow patterns. Int. J. Multiphase Flow 1999, 25, 377–394. (26) Wolk, G.; Dreyer, M.; Rath, H. J. Flow patterns in small diameter vertical non-circular channels. Int. J. Multiphase Flow 2000, 26 (6), 1037– 1061, 2000. . (27) Chen, W. L.; Twu, M. C.; Pan, C. Gas-liquid two-phase flow in micro-channels. Int. J. Multiphase Flow 2002, 28, 1235–1247. (28) Chung, P. M.-Y.; Kawaji, M. The effect of channel diameter on adiabatic two-phase flow characteristics in micro channels. Int. J. Multiphase Flow 2004, 30, 735–761. (29) Coleman, J. W.; Garimella, S. Characterization of two phase flow patterns in small diameter round and rectangular tubes. Int. J. Heat Mass Transfer 1999, 42, 2869–2881. (30) Kreutzer, M. T.; van der Eijnden, M. G.; Kapteijn, F.; Moulijn, J. A.; Heiszwolf, J. J. The pressure drop experiment to determine slug lengths in monoliths. Catal. Today 2005, 105 (3-4), 667–672. (31) Woehl, P.; Cerro, R. L. Pressure drop in monolith reactors. Catal. Today 2001, 69, 171–174. (32) Pehlivan, K.; Hassan, I.; Vaillancourt, M. Experimental study on two-phase flow and pressure drop in millimeter-size channels. Appl. Therm. Eng. 2006, 26, 1506–1514. (33) Taylor, G. I. The rear meniscus of a long bubble steadily displacing a Newtonian liquid in a capillary tube. J. Fluid Mech. 1960, 10, 161–165. (34) Marchessault, R. N.; Mason, S. G. Flow of entrapped bubbles through a capillary. Ind. Eng. Chem. 1960, 52, 79–84. (35) Bretherton, F. P. The motion of long bubbles in tubes. J. Fluid Mech. 1961, 10, 166–188. (36) Goldsmith, H. L.; Mason, S. G. The flow of suspensions through tubes, Part II. J. Colloid Sci. 1963, 18, 237–261. (37) Schwart, L. W.; Princen, H. W.; Kiss, A. D. On the motion of bubbles in capillary tubes. J. Fluid Mech. 1986, 172, 259–275. (38) Lebens, P. J. M.; Stork, M. M.; Kapteijn, F.; Sie, S. T.; Moulijn, J. A. Hydrodynamics and mass transfer issues in a countercurrent gas-liquid internally finned monolith reactor. Chem. Eng. Sci. 1999, 54, 2381–2389. (39) Thulasidas, T. C.; Abraham, M. A.; Cerro, R. L. Dispersion during bubble-train flow in capillaries. Chem. Eng. Sci. 1999, 54, 61–76. (40) Heibel, A. K.; Scheenen, T. W. J.; Heiszwolf, J. J.; van As, H.; Kapteijn, F.; Moulijn, J. A. Gas and liquid phase distribution and their effect on reactor performance in the monolith film flow reactor. Chem. Eng. Sci. 2001, 56, 5935–5944. (41) Heiszwolf, J. J.; Kreutzer, M. T.; van der Eijnden, M. G.; Kapteijn, F.; Moulijn, J. A. Gas-liquid mass transfer of aqueous Taylor flow in monoliths. Catal. Today 2001, 69, 51–55.

(42) Oliver, D. R.; Wright, S. J. Pressure drop and heat transfer in gasliquid slug flow in horizontal tubes. Br. Chem. Eng. 1964, 9, 590–596. (43) Horvath, C.; Solomon, B. A.; Engasser, H.-M. Measurement of radial transport in slug flow using enzyme tubes. Ind. Eng. Chem. Fundam. 1973, 12, 431–439. (44) Bercˇic, G.; Pintar, A. The role of gas bubbles and liquid slug lengths on mass transport in the Taylor flow through capillaries. Chem. Eng. Sci. 1997, 52, 3709–3719. (45) Unverdi, S. O.; Tryggvason, G. A Front-Tracking Method for Viscous, incompressible, Multi-Fluid Flows. J. Comput. Phys. 1992, 100, 25. (46) Tryggvasson, G.; Bunner, B.; Esmaeeli, A. A front tracking method for the computations of multiphase flow. J. Comput. Phys. 2001, 169, 708– 759. (47) Sussman, M.; Smereka, P.; Osher, S. A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow. J. Comput. Phys. 1994, 114, 146. (48) Hirt, C. W.; Nichols, B. D. Volume of Fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 1981, 39, 201–225. (49) Brackbill, J. U.; Kothe, D. B.; Zemach, C. A Continuum Method for Modeling Surface Tension. J. Comput. Phys. 1992, 100, 335. (50) Rider, W. J.; Kothe, D. B. Reconstructing volume tracking. J. Comput. Phys. 1998, 141, 112–152. (51) He, X.; Luo, L.-S. Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation. Phys. ReV. E 1997, 56, 6811–6817. (52) Sankaranarayanan, K.; Shan, X.; Kevrekidis, I. G.; Sundaresan, S. Analysis of drag and virtual mass forces in bubbly suspensions using an implicit formulation of the lattice Boltzmann method. J. Fluid Mech. 2002, 452, 61–96. (53) Taha, T.; Cui, Z. F. Hydrodynamics of slug flow inside capillaries. Chem. Eng. Sci. 2004, 59, 1181–1190. (54) Taha, T.; Cui, Z. F. CFD modelling of slug flow inside square capillaries. Chem. Eng. Sci. 2006, 61, 665–675. (55) Taha, T.; Cui, Z. F. CFD modelling of slug flow in vertical tubes. Chem. Eng. Sci. 2006, 61, 676–687. (56) van Baten, J. M.; Krishna, R. CFD simulations of mass transfer from Taylor bubbles rising in circular capillaries. Chem. Eng. Sci. 2004, 59, 2535–2545. (57) van Baten, J. M.; Krishna, R. CFD simulations of wall mass transfer for Taylor bubbles in circular capillaries. Chem. Eng. Sci. 2005, 60, 1117– 1126. (58) Yang, Z. L.; Palm, B.; Sehgal, B. R. Numerical Simulation of bubbly two-phase flow in a narrow channel. Int. J. Heat Mass Transfer 2002, 45, 631–639. (59) Qian, D.; Lawal, A. Numerical study on gas and liquid slugs for Taylor flow in a T-junction micro channel. Chem. Eng. Sci. 2006, 61, 7609– 7625. (60) Yu, Z.; Hemminger, O.; Fan, L. S. Experiment and lattice Boltzmann simulation of two-phase gas-liquid flows in micro-channels. Chem. Eng. Sci. 2007, 62, 7172–7183. (61) Youngs, D. L. Time-dependent multi-lateral flow with large fluid distortion. In Numerical Methods for Fluid Dynamics; Morton, K., Baines, M., Eds.; Academic Press: New York, 1982; pp 273-285. (62) Renardy, M.; Renardy, Y.; Li, J. Numerical simulation of moving contact line problems using a volume-of-fluid method. J. Comput. Phys. 2001, 171, 243–263. (63) Spelt, P. A level set approach for simulations of flows with multiple moving contact lines with hysteresis. J. Comput. Phys. 2005, 207, 389– 404. (64) Gerlach, D.; Alleborn, N.; Buwa, V. V.; Durst, F. Numerical simulation of periodic bubble formation at submerged orifices under constant inflow conditions. Chem. Eng. Sci. 2007, 2109–2125. (65) Gerlach, D.; Biswas, G.; Durst, F.; Kolobaric, V. Quasi-static bubble formation on submerged orifices. Int. J. Heat Mass Transfer 2005, 48, 425– 438. (66) Buwa, V. V.; Gerlach, D.; Durst, F.; Schlu¨cker, E. Numerical simulations of bubble formation on submerged orifices: Period-1 and period-2 bubbling regimes. Chem. Eng. Sci. 2007, 62, 7119–7132.

ReceiVed for reView May 20, 2008 ReVised manuscript receiVed October 16, 2008 Accepted October 17, 2008 IE800806F