Experimental and Numerical Studies on Microbubble Generation and

Article pubs.acs.org/IECR

Experimental and Numerical Studies on Microbubble Generation and Flow Behavior in a Microchannel with Double Flow Junctions Dongren Liu, Xiang Ling,* and Hao Peng* Jiangsu Key Laboratory of Process Enhancement and New Energy Equipment Technology, School of Mechanical and Power Engineering, Nanjing Tech University, No. 30 Pu Zhu South Road, Nanjing 211816, P. R. China ABSTRACT: Experimental visualizations and numerical simulations were carried out to investigate the bubble generation mechanisms in a microchannel with double flow junctions. Two kinds of flow patterns could be observed from the experiments and simulations: slug-bubble flow and split flow. When 0.00554 < Ca2 < 0.025 and 0.25 < QG/QL < 6, slug-bubble flow can be observed. When Ca2 < 0.00554 and 6 ≤ QG/QL ≤ 8, split flow will occur. Bubble length and velocities were also analyzed experimentally and theoretically. The relationship between bubble length, Ca and flow rate ratio QG/QL was obtained by regression, and the correlation could predict the data within an RMS error of 20%. Bubble velocity data correlated with the drift flux model showed that the distribution parameter of the two flow junctions C01 andC02 had values of 1.0654 and 1.1123, respectively.

1. INTRODUCTION There has been considerable development in microfluidic technology during the past few decades due to its extensive applications in chemical engineering and many other fields. With its high efficiency for heat and mass transfer, microfluidic technology has excellent properties for application in devices such as heat exchangers and mixers. Motivated by these advantages, attention has been paid to the characteristics of fluid flow in microchannels with the goal of amplifying the benefits available from the microscale. In recent years, to meet the needs of chemical and biological production, the growing interest in microfluidics has mainly concentrated on multiphase flow in microchannels with multiflow junctions other than conventional T-junction microchannels (Zhang et al.,1 Bergoglio et al.,2 Lin,3 Wang et al.4). In the study of Zhang et al.,5 two-phase flow in two parallel channels with branches was investigated experimentally with flow patterns and bubble formation mechanisms were observed from the images of a high-speed camera. Relationships between bubble length, Weber number, and the flow rate ratio were obtained by regressing the experimental data. Barreto et al.6,7 studied the pressure drop and flow patterns of two-phase flow in microdevices with seven parallel microchannels, and they found that when the total mass velocity varies from 50 to 300 kg/(m2 s), the flow pattern in each microchannel is different, although the channel configuration is approximately homogeneous. The pressure drop was also measured when the superficial velocities were varied from 0.1 to 3.5 m/s (water) and 0.1 to 34.8 m/s (air). The results showed that both the Reynolds number and void fraction contributed to the pressure drop prediction. In the experiment of Bai et al.,8 a new structural microchannel device was designed. A partially hydrophilic and hydrophobic © XXXX American Chemical Society

glass microchip was used to generate water-in-oil-in-water double emulsions. Wu et al.9 reported a kind of two-phase flow in a microchip with many heart-shaped cells. The results from experiment and CFD simulation showed that momentum exchange was affected by bubble size when the flow rate was lower than 15 mL/min. All the publications mentioned above demonstrated that the novel structure with multiflow junctions had better heat and mass transfer than T-junction or flow focusing junctions, a characteristic that is highly desirable for applications in pharmacy and biology. The flow pattern and the mechanism of bubble formation, especially the characteristics of fluid flow in a microchannel, is still the focus of much research. Many experimental and theoretical analyses were carried out to understand the mechanism of fluid flow and bubble formation at the microscale (Li et al.,10 Yan et al.,11 Santos and Kawaji,12 Cherlo et al.,13 Baround et al.14). Numerous models and correlations were built to summarize and predict the bubble length and velocity (Kawaji et al.,15 Bashir et al.,16 Yamamoto and Ogata17,18). All these theories provide us with a foundation for further study. During the investigation of fluid flow in complicated structures, the flow of a double emulsion in microchannels with double T-junctions and microdevices with double distribution channels were studied as well. The dynamics of bubble breakup in microchannels with double distribution junctions and a divergence were studied experimentally by Fu et al.19 Three-phase flow in microchannels with double T-junctions has been investigated theoretically and Received: Revised: Accepted: Published: A

July 16, 2016 November 5, 2016 November 8, 2016 November 8, 2016 DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

rate to keep a specified flow rate ratio. Gas and liquid superficial velocities were controlled by adjusting the fluid flow rate. When the two-phase flow left the microchannel outlet, it was collected and separated by a water tank. The image of the bubble formation and flow patterns in the microchannel were captured by a highspeed video camera (Phantom, LCR320S, VRI., USA) and a microscope lens with 10× magnification. A high brightness cold light source was used as background illumination. The microchannel with two flow junctions and a rectangular cross section used in the experiment was produced on a poly(methyl methacrylate) (PMMA) sheet. The microchannels were first manufactured by a computer numerical control (CNC) machine on a PMMA sheet. The deviation of the geometry caused by CNC machining has been controlled within ±5%; the surface roughness, Ra, of the microchannel is 10 μm. The PMMA sheet with the microchannel was sealed by another smooth PMMA sheet using the thermal bonding method. The PMMA sheets were put into the vacuum thermal bonding machine and the temperature was increased to the critical temperature to allow the two sheets to bond. This process also increases the surface wettability. The width and depth of all channels were 300 μm. The length of the main channel is 50 mm. Air and water inlets with a width and a depth of 300 μm were distributed symmetrically on both sides of the main channel where the four inlets were normal to the main channel. Each of the inlets was 5 mm long. A schematic diagram of the microchannel and the coordinate system is presented in Figure 2. The microchip

experimentally by Rajesh and Buwa,20 in which the flow regime map was obtained and the air slug lengths were found to be related to the capillary number, Ca, and the Weber number, We. In the work of Wang et al.,21 bubble size laws and flow patterns of three-phase flow in a microchannel with double T-junction were analyzed experimentally and theoretically. The results show that the flow rate ratio is a key parameter in bubble volume prediction. The analyses of the bubble formation mechanism of three-phase flow provide us with a preliminary understanding in this field, with which people are less familiar. In this paper, experimental visualization and numerical simulations were conducted on a microchannel with double flow junctions. This structure was first designed as a precooler for volatile organic compounds (VOCs) recovery. VOCs gas (discontinuous phase) will be cooled by another liquid in the microchannel. As has been reported by Fukagata et al.,22 Dai et al.,23 and Leung et al.,24 heat transfer can be enhanced by slug flow in a microchannel, and the average Nusselt number can reach more than 7 times the Nusselt number for single phase flow in a microchannel. However, the traditional T-junction microchannel cannot meet the requirement of high flow rate in the VOC recovery process. Thus, a microchannel with more inlets was evaluated to solve the problem. Since the bubble length and flow pattern has a significant effect on heat transfer performance, it is very important for us to understand the characteristics of fluid flow in the microchannel. The bubble creation mechanism and the characteristics of two-phase flow in the microchip were also investigated. Furthermore, this structure can also be used for mass transfer in chemical engineering, biological engineering, etc. To facilitate the investigation, three or four phase flow was not considered since two-phase flow was already complicated.

2. EXPERIMENTAL SETUP The schematic diagram of the experimental setup is illustrated in Figure 1. Four syringe pumps were used as stable and continuous sources of gas and liquid flow. The syringe pumps were connected to the microchip by polytetrafluoroethylene (PTFE) tubes. Two streams of gas flow into the microchannel at the same flow rate through two inlets at one side of the main channel and liquid flows into the channel from the other side in the same manner. The liquid flow rates varied with the changes of gas flow

Figure 2. Schematic of the microchannel.

was fixed horizontally on the experimental table to avoid vibration. Bubbles were generated in the microchannel, and then the gas and liquid mixture left the channel from the outlet. Liquid was reserved and the gas was separated in a phase separator. The bubble/slug lengths were measured from the images from the high-speed camera by using image processing software PCC2.6 (Vision Research, Inc. USA). The output resolution was 1024 ×

Figure 1. Schematic diagram of the experimental setup. B

DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

model is used to track the interface between the gas and liquid. The 3D geometry model is established by GAMBIT in which computational domains are discretized by hexahedral cells. Since the 3D model is computationally expensive, a truncated geometry model for the microchannel is built. The width and depth of the channel and four inlets are all W = 300 μm. The length of the gas and liquid inlets are set to 3 W, and the length of the mixing channel is set to 20 W to reduce the computational time. Only several full simulations are conducted for the extreme conditions. 3.2. Governing Equations. The governing equations consist of the incompressible Navier−Stocks (N−S) equation and the continuity equation, coupled with the volume fraction equation. A source term F⃗ correlated with interfacial tension is added to the N−S equation.

1024 with a sample rate of 2700 fps and an exposure time of 50 μs. Water and air were chosen as the working fluids. The physical properties of the working fluids are listed in Table 1. The gas and Table 1. Physical Properties of the Working Fluids liquids

density (kg/m3)

viscosity (Pa·s)

surface tension(N/m)

air water

1.2 998

0.000018 0.0012

0.072

liquid flow rate varied from 0.54 mL/min to 2.16 mL/min and 0.27 mL/min to 6.56 mL/min, respectively, which corresponded to gas and liquid superficial velocities of 0.1 to 0.4 m/s and 0.05 to 1.2 m/s, respectively. The different operating conditions used in this study are listed in Table 2. Each group of experiments was repeated at least five times to ensure the reliability of the experimental data.

∂ρ + ∇ (ρ u ⃗ ) = 0 ∂t

∂(ρu ⃗) + ∇(ρuu⃗ ⃗) ∂t

Table 2. Experimental Conditions QG (mL/min)

QL (mL/min)

QG/QL

UG (m/s)

Ca1

0.54 1.08 1.62 2.16

0.07−2.16 0.27−4.32 0.20−6.56 0.27−8.64

0.25−8 0.25−8 0.25−8 0.25−8

0.1 0.2 0.3 0.4

0.00156−0.00694 0.00313−0.01390 0.00469−0.02100 0.00625−0.02780

= −∇p + ∇[μ(∇u ⃗ + ∇u ⃑ T)] + ρg ⃗ + F ⃗

Ca1 = μL (Q L1 + Q G1)/W 2σ

(1)

Ca 2 = μL (Q L1 + Q G1 + Q L2 + Q G2)/W 2σ

(2)

(5)

The source term F⃗ could be computed from the following: F ⃗ = 2σ

ρ(∇n)̂ ∇αG ρG + ρL

(6)

where α is the volume fraction of the fluid; the subscript L or G refers to gas or liquid; and αG and αL are the volume fraction of gas and liquid, respectively, in each cell. In liquid−gas two-phase flow, when a cell is filled with liquid, the liquid volume fraction, αL, is 1 and the corresponding gas volume fraction, αG, is 0. If the cell is partially filled with liquid and the rest filled with gas, the sum of the volume fraction for all phases is 1:

The capillary number at junction 1 (J1) which is defined as Ca1 and the capillary number at junction 2 (J2) which is defined as Ca2 are expressed as

∑ αi = 1

QL1 and QG1 refers to the gas and liquid volumetric flux injected in to the channel from the first gas and liquid inlets, respectively, and QL2 and QG2 refer to the gas and liquid volumetric flux of the gas and liquid at the second inlets, respectively. W indicates the width of the channel, μ is the viscosity of the liquid, and σ is the surface tension of the fluid. In our experiments, QL1 = QL2, QG1 = QG2, so Ca2 = 2Ca1 and the flow rate ratio QG1/QL1 = QG2/QL2 = (QG1 + QG2)/(QL1 + QL2), so the flow rate ratio in this work will be expressed as QG/QL. The uncertainties of the experiment are dependent on the experimental conditions and the measurement instruments. In this work, the uncertainty analysis is performed using the estimation method. Given a variable R, which is a function of n independent variables Xi = {x1, x2, ..., xn}, the relative uncertainty of R can be expressed as 1/2 ⎡ n ⎛ ⎞2 ⎤ X Δ ΔR i = ⎢∑ ⎜ ⎟⎥ ⎢⎣ i = 1 ⎝ Xi ⎠ ⎥⎦ R

(4)

(7)

n̂ is the interface normal defined by n ̂ = ∇αL /|∇αL|

(8)

The movement of the phase interface in the microchannel is largely affected by the contact angle at the channel wall. To model this effect, the interface normal, n̂, in the cells adjacent to the wall is modified to n ̂ = n ̂w cos θ +

1 − n ̂w 2 sin θ

(9)

where n ̂w is the unit vector normal to the wall and θ is the contact angle. In this study, the contact angle has been set to 36°, which has been used by researchers such as Santos and Kawaji12 and which shows good agreement with the experiment. In a computational cell, the density and the viscosity of the twophase mixture are defined as ρ = αGρG + αLρL (10)

(3)

μ = αLμL + αGμG

where ΔXi are the absolute uncertainties of Xi. By using the above method, the relative errors for the bubble length and the bubble velocity are 2% and 1.4%, respectively.

(11)

The pressure implicit with the splitting of operators (PISO) algorithm is selected for pressure−velocity coupling, the pressure staggering option scheme (PRESTO!) is used for pressure interpolation, and the momentum equation is discretized using the second-order upwind scheme. The interface interpolation is carried out using the Geo-reconstruct scheme. The time step varied from 10−6 s to 10−8 s to maintain convergence. The

3. SIMULATIONS 3.1. Model Geometry. To facilitate the analyses, a series of 3D simulations are conducted using the computational fluid dynamics(CFD) software FLUENT. The volume of fluid (VOF) C

DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 3. Grid independence: (a) effect of mesh number on bubble length; (b) volume fraction of gas and liquid (red is gas and liquid is blue).

effectiveness of the numerical method used in this study has been proven by Kawaji et al.15 and Yamamoto and Ogata.17,18 3.3. Grid Independence and Validation. To examine the grid independence of the simulations, the mesh number is converted from 9500 to 1188000. The flow rates are fixed at QG1 = QG2 = 1.62 mL/min and QL1 = QL2 = 1.62 mL/min, and slugbubble flow can be obtained from the CFD result under this operating condition. The dimensionless bubble length of long slugs created by J1 and J2are chosen for the analyses. The effect of mesh number, N, on bubble length can be seen from Figure 3a. When the mesh number increases from 9500 to 148 500, the dimensionless bubble length increases with increasing mesh number, but when the mesh number varies from148 500 to 1 188 000, the slug length becomes stable and the interface of gas and liquid become clear (shown in Figure 3b). This result indicates that the mesh number will not affect the bubble length significantly when it is greater than 148 500. However, increasing the mesh number will lead to a dramatic increase in computational time. To save computational time but to provide an acceptable accuracy, a mesh number of 352 00 was used for the simulations. The numerical results were validated with the data from the flow visualization experiment. When the flow rate ratio, QG/QL, varied from 0.25 to 6, slug-bubble flow can be obtained both from the experiment and the simulation. Comparison between bubble lengths of slug-bubble flow from experiments and simulations are illustrated in Figure 4, where the solid line represents the dimensionless bubble length from the experiment and the solid dots represent the simulation results. The simulation results are in good agreement with the experimental data, validating the algorithm. The RMS error can be controlled within 10%. The VOF model performed very well on tracking the interface of bubble flow in the microchannels.

Figure 4. CFD validation.

4. RESULTS AND DISCUSSION 4.1. Bubble Formation and Flow Pattern. In the present work, as shown in Figure 2, the first flow junction is depicted as J1, and the second flow junction is defined as J2. To give an overall view of our work, the flow patterns of J1 are also reported in this section. Figure 5 panels a−c are the flow patterns and bubble generation process from the experiment and simulation at different times. The images were taken with a high-speed camera after a steady flow was achieved. A typical bubble generation process in J1 is shown in Figure 5a. According to the results from the experiment and the simulation, slug and round-shaped bubbles were produced when the gas and liquid flow rate was in the ranges of QG1 = 0.54 mL/min to 2.16 mL/min and QL1 = 0.067 mL/min to 8.64 mL/min, respectively. When the capillary number Ca1 is lower than 0.0125, elongated bubbles can be observed, the entire channel is blocked by the growing droplets, and bubble creation is dominated by a squeezing mechanism. As shown in Figure 5c, the bubble formation process will be affected by a shearing mechanism (Chen et al.25) for larger values ofCa1. D

DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Ca2 is lower than 0.014 will the small bubbles created under the squeezing regime appear. It can be seen from the experiment that when the capillary number Ca2 is higher than 0.014, the forming satellite bubble cannot block the channel, liquid film exists during the whole process of bubble formation and the shear stress will affect the whole process of bubble formation. However, when Ca2 is lower than 0.014, the forming satellite bubble will block the entire channel before it detaches from the junction 2 and the squeezing force arising from the liquid from J1 and J2 will act on the satellite bubble tip until it detaches from J2. Therefore, the slug-bubble flow was formed under two flow regimes: squeezing and dripping, which indicates a complete generation cycle including several periods. In the first stage the bubble B1 merged with the growing bubble. Then, the larger bubble B′1, separated from J2, moves forward toward the downstream channel under the squeezing and cutting effect of the continuous phase. At last, one or several small round bubbles are produced in J2 and formed into a satellite bubble of B′1. Another flow pattern is “split flow”. As shown in Figure 7a−e, this flow pattern will be generated when the flow rate ratio QG/ QL is higher than 6 and the capillary number Ca2is lower than 0.00554. In such a flow pattern, the water flow rate is much lower than the gas flow rate, thus the flow regime in both flow junctions is dominated by the squeezing regime when the viscous force is much lower than the surface tension. The main channel will be entirely blocked by the phase injected from inlet G2. When the bubble generated in J1 is passing through J2, it will first merge with the gas slug produced in J2. Then, under the pressure of the continuous phase from L2, the long slug will be split into two smaller bubbles as shown in Figure 7. An arch-type chamber filled with water was formed at the side of the channel. With the increase of volume fraction of water in the main channel, the bubble was squeezed by water in a perpendicular direction. The breakup point is always downstream of the junction 2 due to the movement of the bubble and water from L2. Since the channel has been blocked by the long bubble, the water injected from L2 will be confined by the gas phase and change its flow direction quickly. However, this “split” only occurs once during the bubble generation cycle, which is a little different from the situation described in a previous publication (Wang et al.21). In the study of Wang et al.,21 the breakup point appears at the downstream channel only when the flow is dominated by the jetting regime, and in other situations, the breakup point is at the flow junction. Presumably that resulted from the confinement of the disperse phase when the water was forced to change its flow direction after entering the main channel under the interfacial force. With the continuous feeding of water from L2, the squeezing force increased and the slug was cut into two parts. 4.2. Bubble Shape. Bubble shape has been analyzed by many researchers such as Taha and Cui27 and Fries et al.28 In their studies, the bubbles have a spherical end when Ca is low. However, as Ca increases, the indentation will be seen at the end of the bubble instead of a spherical tail. The cross section of the bubble is nearly square when Ca is lower than 10−3.On the other hand, it will turn to a circle when Ca is higher than 10−1. In addition to Ca, the bubble shape is also affected by the contact angle, which is not involved in our research. The bubble cross section from the numerical results is shown in Figure 8(a−c). The slug cross-sectional area increased with decreasing capillary number Ca2 and decreased with decreasing flow rate ratio, but the corner liquid film can be observed when Ca2 is higher than 0.008. In the work of Fries et al.,28 however, this kind of liquid film will appear when Ca2 is approximately 0.001. The bubble

Figure 5. Bubble generation and flow pattern in J1: (a) typical bubble generation process of experiment and simulations when QG1 = 1.08 mL/ min, QL1 = 1.08 mL/min; (b) experiment and simulation of the flow pattern under shearing regime when QG1 = 1.62 mL/min, QL1 = 3.24 mL/min; (c) experimental and simulation results when QG1 = 1.62 mL/ min, QL1 = 0.202 mL/min.

Thus, most of the experimental conditions were under the regime where squeezing was dominant. The bubble generation can also be divided into several periods. The bubble cap first enters the main channel and then blocks the main channel and grows bigger. Finally, the neck of the bubble splits from J1 and a separated bubble will be observed (Garstecki et al.26) At J2, the bubble generation mechanism is more complicated. The first and the most important kind of flow pattern is “slugbubble flow” (Figure 6f,g). When Ca2 is lower than 0.025 and higher than 0.00554, this flow pattern can be observed with the high-speed camera. As shown in Figure 6f,g, a stable two-phase flow consisting of one gas slug followed by one or several small round bubbles was generated at the downstream channel of the junction 2. Figure 6a−e is a typical formation process of slugbubble flow at J2. In this flow pattern, when a gas slug, B1, created in J1 moves downstream, it will merge with a growing bubble created in J2. During this process, the continuous phase from L2 will squeeze the bubble wall in a perpendicular way, but this force is insufficient to split the gas slug into two small bubbles, so an arch-type curve can be observed in Figure 6c. Another force affecting the process is the squeezing force of the continuous phase from the upstream channel acting on the ends of the gas slug. Under the influence of these forces, B′1 is separated from J2 and becomes an independent bubble. Since the bubble generation frequency in J2 is much higher than that in J1, a round-shaped bubble, B2, will be separated from J2 and became the satellite bubble of B1 before the next gas slug generated in J1 arrives at J′2. With the increase of water flow rate, more than one satellite bubble will be created during one bubble formation period. Most of the satellite bubbles were generated under the effect of shearing force (as shown in Figure 6a−e). Only when E

DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 6. Flow pattern and formation process of slug-bubble flow. (a−e) The bubble formation process of slug-bubble flow (QG1 = QG2 = 1.62 mL/min; QL1 = QL2 = 1.62 mL/min). (f) QG1 = QG2 = 1.62 mL/min. (g) QG1 = QG2 = 2.16 mL/min; QL1 = QL2 = 2.16 mL/min.

⎛Q ⎞ Lb = 1 + ε⎜⎜ G ⎟⎟ W ⎝ QL ⎠

shape observed from the experiment is presented in Figure 8d, and good agreement with the simulation result can be seen in Figure 8e, where the difference between the radius of curvature of the bubble cap and the bubble tail can be observed. In our research, bubbles with a flat tail can be achieved when Ca2 is greater than 0.025, which is much lower than that of Taha and Cui27 and other authors. This is because the bubbles created in J1 will be regenerated in J2 in the microchannel with four parallel inlets. The forces acting on these bubbles during the process of regeneration in J2 is much higher than in J1. Except for the forces from the water which flow into the channel from junction 2, the forces from the liquid upstream must be accounted for as well. These forces lead to a higher interface pressure which affects the bubble shape. 4.3. Bubble Size. Heat and mass transfer performance in microchannels is remarkably affected by bubble and liquid slug length, which has been the focus of many groups. In the work of Garstecki et al.,26 the dimensionless slug length was treated as a function of the flow rate ratio and the fitting parameter ε, expressed as

(12)

where Lb is the bubble length and ε can be defined as

ε=

d W

(13)

where d is the characteristic width and W is the width of the channel. Thus, the fitting parameter ε can be obtained from the experiment. In the research of Garstecki et al.,26 ε = 1, but in Shao et al.,29 they found ε = 0.57 is more suitable for their experimental data. After that, on the basis of this model, a similar correlation with two fitting parameters was proposed by van Steijin et al.:30 ⎛ Q disp ⎞ Lb ⎟⎟ = ε1 + ε2⎜⎜ W ⎝ Q cont ⎠

(14)

Because of its higher accuracy, this model has also been adopted by people engaged in the research of bubble size. F

DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 7. Split flow from visualized experiment and numerical simulations of split flow. (QG1 = QG2 = 1.62 mL/min; QL1 = QL2 = 0.27 mL/min).

Figure 8. Cross section of the gas slug and bubble shape in experiment and simulation. (a−c) Cross section of gas slug(air is red and water is blue). (d) Bubbles with a flat tail when QG2 = 1.62 mL/min and QL2 = 3.14 mL/min. (e) Simulation results when QG2 = 1.62 mL/min and QL2 = 3.14 mL/min.

Figure 10. Comparison between calculated value and experimental data.

Figure 9. Effect of flow rate ratio on bubble length at J1..

et al.26 or Steijin et al.30 As seen from Figure 9, the regression of the experimental data with the linear model yields

In this work, the operating conditions are 0.25 ≤ QG/QL ≤ 6, and 0.00554 ≤ Ca2 ≤ 0.025. Under these operating conditions, slug flow or Taylor flow can be generated at J1. The relationship between dimensionless bubble length L/W and flow rate ratio QG/QL in J1 is illustrated in Figure 9. The dimensionless bubble length increases with increasing flow rate ratio. However, the experimental data did not agree well with the models of Garstecki

⎛Q ⎞ L = 1.06 + 1.70⎜⎜ G ⎟⎟ W ⎝ QL ⎠

(15)

The regressed line agrees well with the data within 20% RMS errors. G

DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

is the difference between the capillary number and viscosity of the fluids. The working fluids used by Steijin30 were air/ethanol, but in our experiment the continuous phase and the disperse phase were water and air. Moreover, the model established by Garstecki et al.26 and van Steijin et al.30 was particularly designed for CaL < 10−2, especially when CaL is in the range of 8 × 10−5 to 8 × 10−3. Equation 12 agrees well with the experiment. CaL is defined as CaL = μLUL/σ, where μ is viscosity, U is the mean speed of the carrier fluid, and σ is the surface tension. In our study, CaL varies from 10−4 to 10−2 where the range of capillary number is wider than in Garstecki’s26 work. In Qian and Lawal’s31 work, two-phase flow in a symmetric Tjunction was studied. They concluded that bubble length is mainly determined by phase hold-up and that the bubble length could be expressed as a function of β, Re, and Ca: Lb = 1.637β 0.107(1 − β )−1.05 Re−0.075Ca−0.0687 d

(16)

where β ≈ UG/(UG + UL) and 15 < Re < 1500. Here, Re is defined as Re = djρL/μL, where j = UG + UL, and Ca is defined as Ca = μLj/ σ. Using Qian’s31 correlation, a better agreement was found between the calculated bubble length and our experimental data. The results are shown in Figure 10. When 0.000278 < Ca1 < 0.01 and 0.09 < β < 0.9, a shorter bubble will be created and the deviations could be controlled within 10%. However, when Ca1 is extremely low (lower than 10−3) and β is relatively high (2 < β < 4), the bubble becomes longer and a much bigger deviation between the experimental data and Qian’s31 correlation can be seen in Figure 11. The estimated bubble length using eq 16 was

Figure 11. Velocity field at J2.

Since the bubble length is highly dependent on channel and inlet geometry, the deviation is probably caused by different orientations of the inlets. In the work of Garstecki et al.26 and Steijin et al.30 the gas and liquid inlets were perpendicular to each other, but in our work the inlets are symmetrically distributed between the mixing channels. Another reason for the discrepancy

Figure 12. Pressure distribution. H

DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 13. Bubble length with flow rate ratio and capillary number Ca2.

Figure 14. Bubble velocities for gas slugs created in J1 and J2.

J1 is much higher than J2 (as shown in Figure 12) and compared with a traditional T-junction, a much higher pressure should be provided. The bubble generation mechanism and flow pattern in J2 is more complicated, so it is difficult to get a correlation to fit all flow patterns to predict bubble size, but we still tried to summarize the rules suitable for some certain circumstances. This research is mainly concerned with the bubble length in slugbubble flow. Compared with other flow patterns, “satellite flow”

too high. The reason is probably that under these operating conditions, the bubble formation mechanism is dominated by the squeezing regime, the downstream channel of J2 has been entirely blocked by gas slugs, and the continuous phase from the upstream channel has been obstructed by the gas slug, resulting in bubble formation being affect in J1. To make this question clear, the velocity field and pressure distribution of the fluid flow in J2 has been studied. As shown in Figure 11, although the fluid flow in J1 will not be affected by the fluid flow in J2, the pressure at I

DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research or “slug-bubble” flow is more regular and easier to be manipulated. The dimensionless slug length at different flow rate ratios and liquid capillary numbers Ca2 for the slug-bubble flow are presented in Figure 13. Although the bubble length increases with flow rate ratio (Figure 13a), the bubble length is not linearly dependent on the flow rate ratio. As seen from Figure 13b, the bubble length in the downstream channel of J2 is also affected by the capillary number. The relationship between dimensionless bubble length and the Weber number, defined as WeG2 = ρGUG22D/σ, is shown in Figure 13c. The distribution of the experimental data shows that the Weber number has no obvious effect on the bubble length. Therefore, the bubble length can be treated as a function of the flow rate ratio and capillary number (as shown in Figure 13d, the solid line represents the regression line of eq 17), and the following equation can be obtained: ⎛ Q ⎞0.589 L Ca 2−0.197 = 2.76⎜⎜ G ⎟⎟ W ⎝ QL ⎠

5. CONCLUSION Experimental visualization and numerical simulations were carried out with a rectangular microchannel which has double flow junctions. The fluid dynamics of gas−liquid two-phase flow in this kind of channel was investigated. The analysis of the bubble fabrication process was performed to produce regular, controllable droplets, which may be useful in engineering applications. From our experimental study, two types of flow patterns can be observed via high-speed photography: slug-bubble flow and split flow. Flow patterns are determined by Ca and the flow rate ratio QG/QL. When 0.00554 < Ca2 < 0.025 and 0.25 < QG/QL < 6, slug-bubble flow can be observed, and when Ca2 is lower than 0.00554, the growing bubble will be split under the pressure of the continuous phase injected from L2 and split flow will be formed. In slug-bubble flow, slugs created in J1 will be regenerated when they pass through the second flow junction. Satellite bubbles are created in J2 under squeezing and shearing regimes. Bubbles with hemispheric heads will appear when Ca2 is lower than 0.025, which is quite different from the research of Taha and Cui27 and Friesh et al.28 In their experiments and simulations, bubbles with a flat tail will appear when Ca is approximately 1.35, but in our research, this phenomenon will appear when Ca2 is higher than 0.025. Moreover, capillary number and flow rate ratio are key parameters that play a significant role in determining bubble length, whereas the Weber number is less important. The distribution parameter C0 was correlated with the experimental data using Zuber and Findlay’s32 model. The results show that when the fluid flow was dominated by the squeezing regime, the distribution parameters were C01 = 1.0654 and C02 = 1.1123.

(17)

The correlation could predict the data within an RMS error of 20%. 4.4. Bubble Velocity. As mentioned above, to simplify the analyses, we only focused our research on the bubble velocities of the bubble-slug flow. In the visualized experiment, the bubble velocity was calculated as UB = ΔDf

(18)

where f is the sample rate of the high-speed camera and ΔD is the distance that the bubble cap travels during 1/f. Figure 14 panels a and b show the velocities of bubbles created in J1 and J2 with respect to the total volumetric flux j, where j = UG + UL. The solid line shows the regressed value of C0 based on the well-known drift flux model presented by Zuber and Findlay:32 UB = C0j + UGj

■

*E-mail: [email protected]. *E-mail: [email protected].

(19)

ORCID

where C0 is the distribution parameter and the drift velocity, UGj is set to zero for the horizontal channel. Many publications about bubble velocity were based on this model. In Warnier et al.33 and Luo and Wang’s34 work, C0 was estimated to be approximately 1.20 and 1.21, respectively. In the work of Kawahara et al.,35 C0 is correlated with different dimensionless numbers and the correlation is predicted as C0 = aBo0.19ReL−0.01WeG 0.01

AUTHOR INFORMATION

Corresponding Authors

Dongren Liu: 0000-0003-3813-6838 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the financial support provided by National Natural Science Foundation of China (Grant No. 51576095) and Natural Science Foundation of Jiangsu Province (Grant No. BK20151539) and Major Collegiate Project of Natural Science Foundation of Jiangsu Province (Grant No. 15KJA480001).

(20)

When we substitute the current parameters into eq 20, the calculated value of C0 is much higher than the value observed in the present work. For the bubbles created in the present work, as seen from Figure 14a,b, the dimensionless bubble length had nearly a linear relationship with the total flow rate ratio. The bubble velocities can be calculated with UB1 = 1.0654j1 and UB2 = 1.1123j2, C01 = 1.0654 and C02 = 1.1123, where j1 = UG1 + UL1, j2 = 2j1. The RMS errors can be controlled within 10% and 20%, respectively. From the experiment, the value of the distribution parameter in J2 is higher than that in junction 1, which indicates that the bubbles created in J2 move faster than bubbles created in J2. The first reason for the discrepancy is presumably that gas and liquid injected from J1 will be blocked at J2 when the bubbles are regenerated in J2, which will reduce the velocity of bubbles created in J1. Another reason is that the forces acting on the bubbles during the process of regeneration in J2 are much higher than the forces acting on the bubbles created in J1.

■

J

NOMENCLATURE B = bubble Bo = bond number Ca = capillary number C0 = distribution parameter ΔD = displacement (mm) F = body force (N) G = gas inlet d = characteristic width (m) j = superficial velocity (m/s) J = flow junction L = channel or bubble length (mm)/liquid inlet n = unit normal DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

(16) Bashir, S.; Rees, J. M.; Zimmerman, W. B. Simulations of microfluidic droplet formation using the two-phase level set method. Chem. Eng. Sci. 2011, 66, 4733−4741. (17) Yamamoto, K.; Ogata, S. Drag reduction of slug flows in microchannels by modifying the size of T-junctions. Int. J. Multiphase Flow 2014, 62, 67−72. (18) Yamamoto, K.; Ogata, S. Effects of T-junction size on bubble generation and flow instability for two-phase flows in circular microchannels. Int. J. Multiphase Flow 2013, 49, 24−30. (19) Fu, T.; Ma, Y.; Funfschilling, D.; Li, H. Z. Bubble formation and breakup mechanism in a microfluidic flow-focusing device. Chem. Eng. Sci. 2009, 64, 2392−2400. (20) Rajesh, V. M.; Buwa, V. V. Experimental characterization of gas− liquid−liquid flows in T-junction microchannels. Chem. Eng. J. 2012, s207−208, 832−844. (21) Wang, K.; Qin, K.; Lu, Y. Gas/liquid/liquid three-phase flow patterns and bubble/droplet size laws in a double t-junction microchannel. AIChE J. 2015, 61, 1722−1734. (22) Fukagata, K.; Kasagi, N.; Ua-arayaporn, P.; Himeno, T. Numerical simulation of gas−liquid two-phase flow and convective heat transfer in a micro tube. Int. J. Heat Fluid Flow 2007, 28, 72−82. (23) Dai, Z.; Guo, Z.; Fletcher, D. F.; Haynes, B. S. Taylor flow heat transfer in microchannelsUnification of liquid−liquid and gas−liquid results. Chem. Eng. Sci. 2015, 138, 140−152. (24) Leung, S. S. Y.; Liu, Y.; Fletcher, D. F.; Haynes, B. S. Heat transfer in well-characterised Taylor flow. Chem. Eng. Sci. 2010, 65, 6379−6388. (25) Chen, X.; Glawdel, T.; Cui, N.; Ren, C. L. Model of droplet generation in flow focusing generators operating in the squeezing regime. Microfluid. Nanofluid. 2015, 18, 1341−1353. (26) Garstecki, P.; Fuerstman, M. J.; Stone, H. A.; Whitesides, G. M. Formation of droplets and bubbles in a microfluidic T-junctionscaling and mechanism of break-up. Lab Chip 2006, 6, 437−446. (27) Taha, T.; Cui, Z. F. CFD modelling of slug flowinside square capillaries. Chem. Eng. Sci. 2006, 61, 665−675. (28) Fries, D. M.; Trachsel, F.; Rohr, P. R. V. Segmented gas−liquid flow characterization in rectangular microchannels. Int. J. Multiphase Flow 2008, 34, 1108−1118. (29) Shao, N.; Salman, W.; Gavriilidis, A.; Angeli, P. CFD simulations of the effect of inlet conditions on Taylor flow formation. Int. J. Heat Fluid Flow 2008, 29, 1603−1611. (30) van Steijn, V.; Kreutzer, M. T.; Kleijn, C. R. μ-PIV study of the formation of segmented flow in microfluidic T-junctions. Chem. Eng. Sci. 2007, 62, 7505−7514. (31) Qian, D.; Lawal, A. Numerical study on gas and liquid slugs for Taylor flow in a T-junction microchannel. Chem. Eng. Sci. 2006, 61, 7609−7625. (32) Zuber, N.; Findlay, J. A. Average volumetric concentration in twophase flow system. J. Heat Transfer 1965, 87, 453−468. (33) Warnier, M. J. F.; Rebrov, E. V.; De Croon, M. H. J. M.; Hessel, V.; Schouten, J. C. Gas hold-up and liquid film thickness in Taylor flow in rectangular microchannels. Chem. Eng. J. 2008, 135S, S153−S158. (34) Luo, R.; Wang, L. Liquid velocity distribution in slug flow in a microchannel. Microfluid. Nanofluid. 2012, 12, 581−595. (35) Kawahara, A.; Sadatomi, M.; Nei, K.; Matsuo, H. Experimental study on bubble velocity, void fraction and pressure drop for gas−liquid two-phase flow in a circular microchannel. Int. J. Heat Fluid Flow 2009, 30, 831−841.

p = pressure (Pa) Q = flow rate (kg/s) Re = Reynolds number t = flow time (s) U = velocity (m/s) W = width of the channel (m) We = Weber number Greek Symbols

α = volume fracture/fitting parameter θ = contact angle of the liquid (deg) ρ = density (kg/m3) μ = viscosity (Pa·s) σ = surface tension coefficient (N/m) β = void fraction Subscripts

1 = first flow junction 2 = second flow junction G = gas L = liquid w = wall

■

REFERENCES

(1) Zhang, B.; Luo, Y.; Zhou, B.; Wang, Q.; Millner, P. D. A novel microfluidic mixer-based approach for determining inactivation kinetics of escherichia coli, o157:h7 in chlorine solutions. Food Microbiol. 2015, 49, 152−160. (2) Bergoglio, M.; Mari, D.; Chen, J.; Mohand, H. S. H.; Colin, S.; Barrot, C. Experimental and computational study of gas flow delivered by a rectangular microchannels leak. Measurement. 2015, 73, 551−562. (3) Lin, Y. Numerical characterization of simple three-dimensional chaotic micromixers. Chem. Eng. J. 2015, 277, 303−311. (4) Wang, X.; Zhu, C.; Wu, Y.; Fu, T.; Ma, Y. Dynamics of bubble breakup with partly obstruction in a microfluidic T-junction. Chem. Eng. Sci. 2015, 132, 128−138. (5) Zhang, L.; Peng, D.; Lyu, W.; Xin, F. Uniformity of gas and liquid two phases flowing through two microchannels in parallel. Chem. Eng. J. 2015, 263, 452−460. (6) Barreto, E. X.; Oliveira, J. L. G.; Passos, J. C. Analysis of air−water flow pattern in parallel microchannels: A visualization study. Exp. Therm. Fluid Sci. 2015, 63, 1−8. (7) Barreto, E. X.; Oliveira, J. L. G.; Passos, J. C. Frictional pressure drop and void fraction analysis in air−water two-phase flow in a microchannel. Int. J. Multiphase Flow 2015, 72, 1−10. (8) Bai, Z.; Wang, B.; Chen, H.; Wang, M. Spatial wettability patterning of glass microchips forwater-in-oil-in-water (W/O/W) double emulsion preparation. Sens. Actuators, B 2015, 215, 330−336. (9) Wu, K.; Nappo, V.; Kuhn, S. Hydrodynamic study of single- and two-phase flow in an advanced-flow reactor. Ind. Eng. Chem. Res. 2015, 54, 7554−7564. (10) Li, X. B.; Li, F. C.; Yang, J. C.; Kinoshita, H.; Oishi, M.; Oshima, M. Study on the mechanism of droplet formation in T-junction microchannel. Chem. Eng. Sci. 2012, 69, 340−351. (11) Yan, Y.; Guo, D.; Wen, S. Z. Numerical simulation of junction point pressure during droplet formation in a microfluidic T-junction. Chem. Eng. Sci. 2012, 84, 591−601. (12) Santos, R. M.; Kawaji, M. Numerical modeling and experimental investigation of gas−liquid slug formation in a microchannel T-junction. Int. J. Multiphase Flow 2010, 36, 314−323. (13) Cherlo, S. K. R.; Kariveti, S.; Pushpavanam, S. Experimental and numerical investigations of two-phase (liquid-liquid) flow behavior in rectangular microchannels. Ind. Eng. Chem. Res. 2010, 49, 893−899. (14) Baroud, C. N.; Gallaire, F.; Dangla, R. Dynamics of microfluidic droplets. Lab Chip 2010, 10, 2032−2045. (15) Kawaji, M.; Mori, K.; Bolintineanu, D. The effects of inlet geometry and gas-liquid mixing on two-phase flow in microchannels. J. Fluids Eng. 2009, 131, 041302−1−041302−7. K

DOI: 10.1021/acs.iecr.6b02712 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX