Inertial manipulation of bubbles in rectangular microfluidic channels

Inertial microfluidics is an active field of research that deals with crossflow positioning of the suspended entities in microflows. Until now, the majority of the studies have focused on the behavior of rigid particles in order to provide guidelines for microfluidic applications such as sorting and filtering. Deformable entities such as bubbles and droplets are considered in fewer studies despite their importance in multiphase microflows. In this paper, we show that the trajectory of bubbles flowing in rectangular and square microchannels can be controlled by tuning the balance of forces acting on them. A T-junction geometry is employed to introduce bubbles into a microchannel and analyze their lateral equilibrium position in a range of Reynolds (1 < Re < 40) and capillary numbers (0.1 < Ca < 1). We find that the Reynolds number (Re), the capillary number (Ca), the diameter of the bubble (D[combining macron]), and the aspect ratio of the channel are the influential parameters in this phenomenon. For instance, at high Re, the flow pushes the bubble towards the wall while large Ca or D[combining macron] moves the bubble towards the center. Moreover, in the shallow channels, having aspect ratios higher than one, the bubble moves towards the narrower sidewalls. One important outcome of this study is that the equilibrium position of bubbles in rectangular channels is different from that of solid particles. The experimental observations are in good agreement with the performed numerical simulations and provide insights into the dynamics of bubbles in laminar flows which can be utilized in the design of flow based multiphase flow reactors.


Introduction
Multiphase flows are of critical importance in systems which transcend the boundaries of a single discipline. In thermal engineering, phase change heat sinks provide higher heat transfer rates and uniform surface temperatures compared to their single phase counterparts. 1 A large number of studies have investigated the performance of these heat sinks for electronic cooling. [2][3][4][5] In the field of reactor engineering, variants of reactors including trickle-bed, 6 fluidized bed, 7 and bubble column reactors 8 rely on multiphase flows to enhance mass transport which requires a thorough understanding of their complex physics. Moreover, in most electrochemical sys-tems, gaseous products evolve at the interface of the electrodes and the electrolyte, including important processes such as water electrolysis, 9 chlor-alkali 10 and aluminum production. 11 Better understanding and control of the phenomena associated with multiphase flows will contribute to improving the performance of these devices. For instance, many electrochemical reactors use ion conducting membranes or porous separators between the anodic and cathodic half cells. Recently, we have shown that by tuning the flow rate of a liquid electrolyte in a microfluidic electrolyzer, it is possible to Lab Chip, 2018, 18, 1035-1046 | 1035 This journal is © The Royal Society of Chemistry 2018 a Optics Laboratory, School of Engineering, Swiss Federal Institute of Technology control the trajectory of the generated bubbles in a membrane-less configuration. 12 Removing the solid membrane in these systems decreases the ohmic losses and provides flexibility in the selection of catalysts in addition to a simplified cell design. [13][14][15][16] Controlling two-phase flows by topography induced variations of the bubbles' surface energy is another way of eliminating the need for a functional membrane. 17 In this paper, we implement experimental and numerical tools to investigate the dynamics of bubbles flowing in microfluidic channels. By doing so we are able to determine the key parameters such as the diameter of the bubble, the channel geometry, and the capillary and Reynolds numbers that influence the position of bubbles in the microchannels and provide insights into how to control the bubbles' motion in the flow.
A large body of research in the field of inertial microfluidics has been dedicated to controlling the position of particles in fluidic channels since the discovery of the Segre- image of the bubble formation at the T-junction, and (d) a cross-sectional representation depicting how the radius of the bubble, r, and the bubble distance from the wall, d Y , are measured at 6 mm downstream of the T-junction. (e) Experimental setup: the liquid is injected using a syringe pump and a constant pressure pump is used to introduce the gas phase into the device. Motion of bubbles is recorded with a high speed camera in the reflection mode. Silberberg effect. 18,19 Most of the effort in this field is focused on the control of solid particles with diameters smaller than the size of the channel; typically, by one order of magnitude. [20][21][22][23][24][25][26][27][28][29][30] Although there are similarities between the behavior of rigid particles and bubbles, deformation of the bubbles [31][32][33][34] complicates the understanding of the physical phenomena evidenced in gas/liquid flows. Understanding the active forces affecting bubble motion is necessary in order to explain its behavior. Important forces at moderate Re (1 < Re < 100) are the wall-induced lift force, the shear gradient lift force, and the deformation-induced lift force. These forces are shown schematically in Fig. 1. In this range of Re numbers and length scales, the buoyancy lift force can be neglected since the inertial forces dominate by at least one order of magnitude. The wall-induced lift force originates from the asymmetry in the flow field around the bubble and its direction is towards the center of the channel. 35 This force decreases as the bubble moves further away from the wall since the flow asymmetry decreases. The shear gradient lift force is a result of the non-uniform velocity profile and pushes the bubble towards the wall. 32,36 The non-uniform velocity profile induces unequal relative velocities on the two sides of the bubble and as a result, it changes the pressure distribution around the bubble. In the parabolic velocity profile, the relative velocity of the flow to the bubble is higher at the closer edge to the wall resulting in a smaller pressure at this side. Therefore, the bubble moves towards the wall where the magnitude of the relative velocity is the highest. 37 Deformation of the bubble induces a force towards the centerline where the bubble can minimize its energy by reducing its deviation from the original spherical shape and gain a symmetric shape. [38][39][40] This force is considered as a deformation-induced lift force. The balance of these forces determines the equilibrium bubble position in the crossflow direction.
Numerical and experimental investigations 31,33 have confirmed the presence of the deformation-induced, wall and inertial lift forces on the droplets moving in microchannels. As will be discussed in the numerical simulation section, the    32 In their research, the bubble position is determined by balancing the buoyant force with hydrodynamic lift forces, while at higher Re, the balance of forces acting on the bubble is different. In the study of Rivero-Rodriguez and Scheid, 40 the motion of an unconfined bubble in circular microchannels is investigated analytically and numerically. Although their results confirm the effect of inertial and deformation-induced lift forces, the motion of the bubble in the rectangular microchannels is not included in their study.
In this paper, we investigate the bubble motion at moderate Re (1 < Re < 40). We find that by moving towards larger Re, the inertial lift force starts to dominate and gradually changes the equilibrium position. The bubble position in this range of Re is determined by balancing the inertial lift forces with the wall and deformation-induced lift forces, as the effect of the buoyant force is negligible. We use experiments and numerical simulations to investigate systematically the interplay of the bubble diameter, Re, Ca, and the channel ge-ometry in determining the steady state position of the bubbles. We demonstrate that in the channels with aspect ratios higher than one or at high Re, the bubbles stay close to the wall while increasing the Ca or the bubble diameter causes the bubbles to migrate towards the centerline. We discovered that the behavior of the bubble in the rectangular channels is different from the solid particles and this behavior can be modified by using surfactants.

Experimental setup
A T-junction geometry is fabricated to introduce nitrogen bubbles from a narrow side channel into the wide main channel as shown in Fig. 2. Microfabrication of the channels starts on a double-sided silicon (Si) wafer with a 1.5 μm silicon dioxide (SiO 2 ) layer on both sides. Two photolithography    steps are performed on the topside of the wafer to enable different channel depths for the main and side channels. A third photolithography step is used on the backside to make through holes as fluidic access ports. Channels and ports are fabricated using deep reactive ion etching (DRIE). After etching, a 100 nm layer of silicon dioxide is grown on the silicon wafer. Finally, channels are sealed by bonding the silicon wafer to a Borofloat 33 glass wafer via anodic bonding. Fig. S1 (ESI †) illustrates the complete microfabrication procedure. In order to study the effect of geometry, three sets of channels are fabricated with a depth of 50 μm and widths of 50, 150, and 200 μm. The gas channel is perpendicular to the main channel and its depth and width are both 5 μm. The main channel's length is 16 mm and the gas channel connects to the main channel at the midpoint along the direction of flow. There are two inlets for injecting the liquid and the gas and one outlet for both phases. The flow rate of the liquid and the pressure of the gas are regulated to produce different bubble radii in a wide range of Re values, while Ca is changed by employing different liquids as the carrier phases. Table 1 lists the properties of water, ethanol (Thommen-Furler AG), and isopropanol (Fisher Scientific UK) which are used as carrier phases. The surface tension coefficients of the liquids are measured using a KRUSS DSA-30 drop shape analyzer using the pendant drop method. 41 Gas bubbles are generated at the side of the channel. A syringe pump (Cronus Sigma 1000 Series) and a constant pressure pump (Elveflow OB1 MK3) are used for introducing liquid and gas phases, respectively. Bubble properties are measured at a position 6 mm downstream of the T-junction. Videos are recorded using a high-speed camera (Photron FASTCAM Mini UX100) with a frame rate of up to 10 000 fps. The schematic of the experimental setup is depicted in Fig. 2e. The recorded video is used to measure the bubble's distance to the side and the bubble radius in each frame. Diagrams of the bubble position with the error bars representing the standard deviation from the average value are shown in Fig. S2 (ESI †). Based on the measured values, the following parameters are used to define the bubble position for convenience: (1)  where w, r, d Y , and d Z are the width of the channel, the radius of the bubble, the minimum distance between the bubble's edge to the wall in the y-direction, and that in the z-direction, respectively. According to these definitions, Ȳ and Z are unitless quantities that vary from zero to one as the bubble moves from the wall towards the center. Furthermore, it should be noted that the diameter of the bubble is normalized to the channel width: (3)

Numerical simulation
We carry out simulations in COSMOL Multiphysics. It resolves the fluid interface by the arbitrary Lagrangian Eulerian (ALE) moving mesh technique, where the interface is represented by exact mesh lines body-fitted with the bubble.
The nondimensional incompressible Navier-Stokes equations solved using COMSOL are: (4) where p is the pressure and U is the velocity. Re denotes the ratio of the inertial forces to the viscous forces as: (5) where μ l is the viscosity of the liquid, ρ l is the density of the liquid, U cl is the centerline velocity of the underlying flow, and R is the radius of the bubble at rest. The boundary condition at the interface can be written as: (6) where n is the unit normal vector to the bubble interface, I is the unit tensor, and the l and g subscripts denote the liquid and the gas phases, respectively. The capillary number Ca compares the viscous forces against the capillary forces: where σ is the surface tension at the interface. Because the viscosity ratio , eqn (6) can be simplified as follows: Since the volume of the bubble is supposed to be constant, the pressure p g can be obtained by imposing an additional constraint where the surface integral of the normal velocity on the bubble interface is zero: We perform the simulations in the reference frame that has the same streamwise translational velocity U bubble of the bubble. In this frame, the velocity imposed on the wall is U bubble e x . At the inlet, the analytical velocity profile (shifted by −U bubble ) of a rectangular duct is applied; hence the mean inlet velocity is (U mean − U bubble )e x , where U mean is the average velocity of the underlying flow. The advantage of adopting this frame is that the bubble only translates in the spanwise direction, viz., in the plane normal to the streamwise direction and is hence not advected out of the computational domain but still allowed to freely rotate and deform. This also significantly alleviates the need for expensive remeshing and interpolation which have been the main hurdles of ALEbased simulations. Note that U bubble is a time-dependent unknown obtained at every time step, and it is subject to an additional constraint where the surface integral of the streamwise velocity on the bubble interface is zero. We launch the simulation by initially putting a spherical bubble slightly away from the center of the domain. After a transient period (see Fig. 3), the bubble will reach an equilibrium state, where U bubble , the shape and the lateral position of the bubble do not vary with time. We do not attempt to resolve the transient phase of the dynamics, because in this chosen reference frame, the underlying flow changes the magnitude with time, introducing an artificial acceleration term that will influence the transient dynamics. Luckily, we only focus on the equilibrium configuration of the bubble when it reaches a steady velocity U steady bubble and hence the artificial acceleration term vanishes without influencing the equilibrium dynamics. To confirm that the equilibrium configuration is independent of the transient dynamics, we have conducted for certain cases a new computation where the velocities towards the wall and in the inlet are shifted by a time-invariant magnitude, U steady bubble . The new equilibrium configuration is the same as the previous one, confirming that our strategy is appropriate.
We span a wide range of the three nondimensional numbers Re, Ca, and D for a thorough computational investigation. Fig. 3 shows the computational setup and results. The velocity of the bubble is monitored as a function of time to make sure that the solution converges to the steady state. To verify the conservation of the bubble's volume, the temporal evolution of the bubble's volume is plotted in Fig. 3e. After verifying that the equilibrium state has been reached, the bubble position and the velocity profile around the bubble are extracted and shown in Fig. 3f and g, respectively.

Validation of simulations
Recently, COMSOL Multiphysics has been used for simulation of droplets in circular geometries where the results are in good agreement with the asymptotic analysis and the boundary integral method which shows its capability in resolving the droplet motion. 42 Table S1 † shows the results of the numerical simulations for different grid sizes. The grid parameter sizes of the third simulation are optimal and chosen for the numerical study. In this section, for validation of the numerical simulations with the experiments, simulations are performed over the same range of Re, D, and Ca as in experiments. Experimental data are selected for both of the cases where bubbles go to the center or to the wall. Also, data are chosen from two different types of liquids to probe the impact of Ca. The comparison between the two sets of results is summarized in Table 2 with isopropanol as the carrier fluid and Table 3 with ethanol as the carrier fluid. In the tables, Ȳ and Z are calculated using eqn (1) and (2). At Re ≤ 1.29, the predicted bubble position from the simulations does not correspond to the experimental values. This discrepancy can be explained by considering the possibility that multiple equilibrium positions exist in the channel. 33 As a result, it is possible that the simulation chooses one of the equilibrium positions and the experiment the other. More importantly, in the experiments, a row of bubbles can be seen, while a single bubble motion is simulated numerically. At low Re, the range of influence of a bubble increases; promoting inter-bubble interactions that can affect the lateral position of the bubble. 43 At higher values of Re and Ca, excellent agreement can be seen between simulations and experiments. Table 4 shows the numerical results for the equilibrium position of the bubbles with different initial positions. The results show that the bubble, starting from different positions, reaches the same equilibrium position. However, this result cannot be considered as a general statement since at low Re, a bubble can possibly have multiple equilibrium positions as seen for solid particles. 44

Results
The influence of Re, Ca, the bubble diameter and the channel geometry on the positioning of the bubbles merits a detailed discussion. In the following subsections, each parameter's role is discussed and in the last part, the equilibrium position of the bubbles is compared with the solid particles. This journal is © The Royal Society of Chemistry 2018

Effect of the Reynolds number
The Reynolds number determines the balance of the different forces acting on the bubble. The shear gradient lift force to-wards the wall depends on the relative velocity at the two sides of the bubble. By increasing the Re, the difference in relative velocity at the two sides of the bubble increases and  as a result, the shear gradient lift force increases. The higher relative velocity near the wall generally translates to lower pressure on the side of the bubble facing the wall; therefore, this force dominates at higher Re, pushing the bubbles towards the wall. Fig. 4 shows the effect of the Re in the 50 μm wide channel. At small Re values, the bubble equilibrium position is at the centerline. Fig. 5 shows the numerical results for the bubble's lateral position at different Re values for Ca = 0.1 and D = 0.4. In each simulation, the bubble's initial position is at (Ȳ, Z) = (0.7667, 0.6). As the Re increases, the bubble moves closer to the wall and at Re larger than 20, the bubble moves to the corner of the channel on the diagonal.

Effect of the capillary number
The capillary number signifies the ratio of viscous forces to the surface tension forces. At large Ca, the viscous forces are stronger than the capillary forces and a bubble under this condition undergoes larger deformation. The minimum surface energy for the bubble is achieved in the spherical shape. This condition is best met when the bubble is in the middle of the channel where the velocity profile is symmetric and the velocity is the same on both sides of the bubble. As a result, the flow exerts stronger deformation-induced lift force on the bubble towards the center of the channel. Fig. 6 shows the bubble position at different Ca values. In this figure, the con-tinuous phases in the first and second rows are isopropanol and ethanol, respectively. The Ca is larger in the first row and as a result, the distance between the bubble and the wall is larger in the first row. Numerical results at different Ca values with Re = 10 and D = 0.4 are shown in Fig. 7. The bubble equilibrium position is closer to the centerline at larger Ca. At Ca higher than 0.9, the bubble moves to the centerline.

Effect of the bubble's diameter
Lateral forces acting on the bubble originate from unequal pressure and velocity around the bubble. These parameters change as the diameter of the bubble varies. Fig. 8 presents the effect of bubble diameter on its equilibrium position. It can be observed that the larger bubble is closer to the channel's centerline. Fig. 9 shows the results of the simulations for various bubble diameters at constant Re = 10 and Ca = 0.1. The equilibrium position of a bubble smaller than half of the channel is between the diagonal line and the center of the sidewall. On the one hand, such a bubble moves towards the center of the sidewall similar to a small solid particle. 45 On the other hand, the deformation-induced lift force pushes the bubble towards the diagonal line where the shear stress is smaller as shown in Fig. S5. † Therefore, as the diameter of the bubble decreases, the deformability of the bubble decreases and the equilibrium position of the bubble moves closer to the center of the sidewall. When the bubble diameter becomes larger than half of the width of the channel, the bubble moves along the diagonal and by increasing the diameter of the bubble, eventually, it moves to the centerline of the channel. Larger bubbles feel stronger wall-induced lift force towards the center and besides that, when the diameter of the bubble exceeds half of the channel width, the shear gradient lift force towards the wall decreases. Therefore, as the bubble becomes larger, forces towards the center dominate and the shear gradient lift force towards the wall diminishes resulting in positioning of the bubble at the centerline.

Effect of aspect ratio of the channel
As discussed in the previous section, bubbles with diameters larger than half of the channel width with a square geometry move towards the center. In the channels with aspect ratios larger than one, the bubble diameter can be larger than half of the height and smaller than half of the width of the channel. Therefore, the bubble moves towards the center in the  height direction and moves towards the wall in the width direction. Fig. 10 shows the experimental and numerical results for channels with different aspect ratios. In the square channel, since the bubble diameter is larger than half of the channel's width, it moves to the centerline. In the channels with aspect ratios larger than one, the equilibrium position of the bubble changes from the centerline to the narrower side of the channel. Fig. 11 shows the experimental results for channels with three different aspect ratios. In the square channel, bubbles can be in the center of the channel or close to the wall based on their diameters. As discussed in section 3, for bubbles with larger diameter, the forces towards the center increases and push the bubbles to the center while smaller bubbles move to the wall. On the other hand, in the rectangular channels, the equilibrium position of the bubbles with the same diameters is close to the narrow side of the channel.

The difference between equilibrium positions of bubbles and particles
As shown in Fig. 11 for the considered parameter range, the equilibrium position of the bubble in a rectangular channel is close to the shallower edge, whereas it is known that the equilibrium position of a rigid particle in a rectangular channel is in the center of the channel's longer wall. 35,[45][46][47] The different behavior of the bubble and the solid particle is due to the existence of a rigid particle-fluid interface in the former and a fluid-fluid interface in the latter case. To have a better comparison between the rigid particles and bubbles, it is possible to add a surfactant to the carrier fluid to make the bubbles behave similarly to rigid particles. 48,49 Fig. 12 shows the equilibrium position of particles when the carrier fluid is water and the equilibrium position of bubbles when the continuous phase is water or a solution of water and a surfactant (a 1% volumetric aqueous solution of Triton X-100) in the channel with 150 μm width. This figure demonstrates that even with a small amount of the surfactant, the bubbles move to the center of the wider channel while at the same Re, bubbles in surfactant-free water stay close to the shallower wall. By adding the surfactant to the liquid, molecules of the surfactant are adsorbed at the interface, making it a no-slip surface. This phenomenon reduces the interaction of the two fluids at the interface and makes the bubbles behave similarly to the rigid particles. Fig. 13 shows the schematic of the equilibrium position of clean bubbles, contaminated bubbles, and solid particles in the cross section of the channel. A clean bubble goes towards the shallower wall while a solid particle or a contaminated bubble stabilizes at the center of the wider sidewalls. 45,46 6. Transition of the bubble position

Conclusion
Monodisperse bubbles' motion in rectangular microchannels was investigated in order to observe their crossflow equilibrium position. The diameter of the bubble, the aspect ratio of the channel, and the Reynolds and capillary numbers are found to be the critical parameters. Analysis of the results reveals the pattern of the changes in the bubbles' position. Both numerical and experimental data show that a high Re flow pushes the bubbles towards the wall. Moreover, bubbles in the channel with an aspect ratio higher than unity move to the shallower sidewall and by increasing the aspect ratio of the channel their distance from the center increases. On the other hand, by increasing the capillary number, the bubble moves close to the center and for bubbles with D = 0.4 in a square channel and at capillary numbers higher than 0.9, their equilibrium position is at the centerline. Furthermore, bubbles that have a radius smaller than half of the width of the channel move to the wall while bubbles with a radius larger than half of the channel's width move to the centerline. Finally, the comparison between the equilibrium position of a bubble and a solid particle demonstrates that the shallower wall of the rectangular channel is the equilibrium position of the bubble while the solid particle moves to the center of the wider sidewall. This study provides insights into the effective parameters in the bubble motion and the methods we presented can be employed in designing multiphase microfluidic devices. Future studies on the interactions of bubbles introduced from multiple side channels with various sizes will shed more light on the complex dynamics of bubbly flows.

Conflicts of interest
There are no conflicts of interest to declare.