Modeling the motion of a Taylor bubble in a microchannel through a shear-thinning fluid

Applications of multiphase flows in microchannels as chemical and biological reactors and cooling systems for microelectronic devices typically present liquid slugs alternated with bubbles of elongated shape, the Taylor bubbles. These occupy almost entirely the cross-section of the channel and present a hemispherical front and a liquid layer, the lubrication film, which separates the gas from the tube wall. The Taylor bubble perturbs the surrounding fluids activating many transport mechanisms in the proximity of the gas-liquid interface; therefore, the bubble motion significantly influences the heat and mass transfer rates. Although many works deeply investigate the bubble hydrodynamics in Newtonian fluids, the knowledge about the relation between bubble hydrodynamics and rheological properties is insufficient, and studies where the continuous phase exhibits a shear-thinning behavior are missing. Our numerical analysis tries to fill this gap by investigating the motion of a Taylor bubble in a non-Newtonian shear-thinning fluid, modeled by the Carreau viscosity model. First, we validate the results against the Newtonian case and a recent theory for shear-thinning fluids (Picchi et al., Journal of Fluid Mechanics, 2021, 918). Then, we investigate the bubble hydrodynamics far from the validity range of the current models. Finally, we study the scaling of the bubble velocity and lubrication film thickness, extending the current theory to shear-thinning fluids. ⇥

Microchannels present a rectangular or circular cross-section of a sub-millimeter hydraulic diameter [8], and the viscous and surface tension forces dominate over the inertial and buoyancy forces. Since this regime could produce long di1usion times and low mixing performance [9,10]; multiphase flows can enhance the heat and mass transfer rates [11,12,13]. ⇥ Slug flows are frequent in microfluidics [14,15], and consist of liquid slugs, alternated to bubbles of elongated shape, Taylor bubbles, that occupy almost entirely the cross-section of the channel, and a thin liquid layer, named lubrication film, separates the gas from the tube wall. Since the transport mechanisms change in the di1erent regions of the bubble [16], predicting the bubble shape is crucial to estimate the heat and mass transfer rates. ⇥ First studies on elongated bubble hydrodynamics [17,18] scale the bubble speed with the capillary number (Ca=⌅U/⌥) showing that the bubble moves faster than the mean velocity of the flow by the e1ect of the lubrication film. The theory by Bretherton [19] models the curvature of the bubble front and the lubrication film thickness in relation to Ca.
The most recent works model the bubble profile by more sophisticated techniques to extend Bretherton's theory at broader operating conditions as higher Ca and square-section capillaries [20,21,22]. ⇥ All works mentioned above study the bubble motion in a Newtonian fluid. Shearthinning fluids are non-Newtonian fluids that exhibit a pure viscous behavior [23]; their deformation depends only on the current value of shear stress (the duration of the deformation is irrelevant); the viscosity is maximum at zero shear rate, and it gradually diminishes to the minimum when the shear-rate tends to infinite. These fluids are typical of biological fluids [24,25].⇥ The hydrodynamics of a Taylor bubble in a shear-thinning fluid presents more complexities than in the Newtonian case [26] because the fluid viscosity changes in the bubble surrounding. Although numerous works study the bubble motion in Newtonian fluids, the available literature poorly investigates the relations between the rheological properties and bubble profile, and studies where the continuous phase exhibits a shearthinning behavior are currently missing. ⇥ Our numerical analysis aims to fill this gap. We study the motion of a Taylor bubble in a shear-thinning flow at a broad range of operating conditions, including cases with low and high shear rates. To properly model the bubble motion, we need a constitutive model for the nonn-Newtonian fluid. The power-law is the simplest rheological model, and it finds wide application in microfluidics [28,29,30,31]; however, this model presents poor accuracy when applied to a shear-thinning fluid, especially at a low shear rate. Thus, we solve the Carreau constitutive model [32], which includes an additional term for the viscosity at zero shear rate, and, therefore, it overcomes the limitations of the power-law viscosity model [33]. The solution of the Carreau model introduces further novelty in our works because the application of this constitutive equation in microfluidics is scarce: the current studies on the bubble hydrodynamics in Carreau fluids [34,35,36,37,38,39] limit their results to a single fluid (carboxymethylcellulose) and are valid below the tested operating conditions.⇥ The case study is a planar microchannel. Starting from simulating the bubble motion in a Newtonian fluid, we check the independence of the solution from the grid and then validate the numerical results against experimental data from the literature. After the validation of the computational domain, we simulate the bubble motion in shear-thinning fluids and compare the results with our recent theory [27], valid only for low Ca. Finally, we investigate the bubble hydrodynamics out of the validity range of the theoretical model. Our final goal is to predict the lubrication film thickness and bubble profile in a wide range of Ca.⇥ 2 Materials and method⇥

Numerical domain⇥
The numerical domain is a planar microchannel of rectangular cross-section (Fig.1). The length L of the channel is 500=m, and the radius R is 50=m. As boundary conditions, we impose the symmetry along the horizontal semi-axis of the channel and a velocity profile at the inlet. In the Newtonian cases, the inlet profile is semi-parabolic:⇥ Where Umax and Umean are the center-line and average fluid velocity, respectively. In shearthinning fluids, the profile is obtained from the solution for Carreau fluids presented in [40]. The no-slip condition sets the velocity of the fluid to zero at the downward wall. The outlet presents no constrains.
We solve a structured grid of nx?ny squared elements. Similar works [41,42] show that the lubrication film thickness is particularly sensitive to the number of elements along the vertical direction (ny). Thus, we fix nx=800 in all simulation cases and vary⇥ny⇥to check the grid independence of the solution. The vertical dimension of each element changes along the radius, increasing the number of mesh in the proximity of the tube wall to optimize the calculation of the lubrication film.
At time zero, the channel contains the liquid phase and a bubble of 200=m length, positioned at 20=m from the inlet. The bubble rear and front are hemispherical with a radius of 40=m, and we calculate the initial pressure di1erence between the bubble and the fluid by the Laplace equation [43].

Governing equations
We use the interFoam solver of the OpenFOAM-v8. The solved equations represent a system of two immiscible, incompressible, and isothermal fluids, and the interface tracked by the volume of fluid (VOF) technique. Balance equations of mass and momentum read: Where U is the fluid velocity, p is the pressure, g is the gravitational acceleration, and F ⌥ is the surface tension force. The conservation equations include the contribution of both phases by averaged fluid properties: The term ⇡ defines the ratio between the volume occupied by the primary phase over the total cell volume, and it varies from 0 and 1 at the interface. In the VOF, the volume fraction ⇡ is a passive scalar transported by the flow: The discrete solution of this equation needs an appropriate numerical scheme to limit numerical di1usion and ensure the boundedness of the G between 0 and 1. First order schemes preserve the boundedness but generate high di1usion at the interface, the second order schemes are less di1usive but they generate oscillations of the volume fraction out from the definition interval. In our analysis, we solve the transient term by the first order implicit Euler scheme, and the convection term by the Gaussian integration with a linear interpolation (a second order scheme); the VanLeer limiter [44] guarantees the boundedness of G. Furthermore, the interFoam solver includes the Multidimensional Universal Limiter with Explicit Solution (MULES) algorithm [45], that balances the e1ect of numerical di1usion by an artificial compression term in the momentum equation: The compression term is non-zero only at the interface due to the multiplying factor ⇡(1−⇡), and the compression velocity U r is: where S is the area of the cell face multiplied by the outward vector ns normal to face, and the term "=U·S is the volume flux accross the cell face. The constant C ⌫ is a tuning value which controls the compression; a value of 0 defines no compression, while an higher value ensure a sharp interface. In our simulation, we set C ⌫ =1 because higher values increase the error of the interface curvature and smearing; this choice is common for simulating surface tension dominated flows (see [46,47]).

Operating conditions and dimensionless groups
As first step to investigate the bubble motion in the Carreau fluid, we compare our numerical results with the recent lubrication model for shear-thinning fluids in [27]. In this theory, the liquid phase follows the Ellis constitutive equation [33], and we study the bubble motion at di1erent degrees of shear-thinning, measured by the Ellis number El. As defined in [27], El is the ratio between the characteristic shear-rate of the fluid and the representative shear-rate in the lubricating film. When El tends to infinite, the viscosity model reduces to Newtonian; at small El, the shear-thinning e1ects dominate.
The above theory is valid as Ca tends to zero. Thus, we use it to validate our results at low capillary numbers, and then we investigate broader operating conditions by pure numerical solutions of the bubble motion. We simulate the Taylor bubble in a shearthinning fluid at three orders of magnitude of Ca. At each capillary number, we investigate four rheologies (i.e., four combinations of the Carreau's input parameters), each one corresponding to a specific El value, so to space from Newtonian to pure shear-thinning behavior. See the simulated rheology curves in Fig.2 and the input parameters of the Carreau model in Table1.

Grid check and validation of the numerical domain
This section presents the grid independence check and validation of our numerical solutions against current theories for Newtonian fluids.
We check the grid independence by the film thickness calculated at Ca= 10 C3 by grids with di1erent n y values. Simulations at such Ca demand the highest level of accuracy (and computational e1orts) because the lubrication film typically is two orders of magnitude thinner than the tube radius when Ca tends to zero [19,41]. Results (Fig.3) show the progressive stabilization of H∞ in finest grids (n y =225C250), where it converges to the theoretical value from Bretherton [19], with a percentage error within 6-8%. The accuracy augments by the number of mesh; however, too high n y values lead to unfeasible computational costs. Therefore, we locally refine the grid near the tube wall by di1erent expansion rates, namely, the ratio of the coarsest over the thinnest mesh length along the radial direction (' 1 /' ny in Fig.4). Results show the most severe expansion rates (' 1 /' ny N105) resulting in the lowest percentage error (5.8%). The validation compares the dimensionless film thickness obtained from simulations against the theory by Bretherton and results show an excellent agreement at low capillary numbers (Fig.5), with a maximum di1erence of 8.7%. Our model diverges from theory when increasing the Ca, and it presents a maximum di1erence of 20% at Ca=2.4?10 C2 . This trend is reasonable because Bretherton's theory, valid when Ca tends to zero, loses accuracy when the capillary number increases. At high Ca, we compare our results against the correlations proposed by Balestra et al. [48]. The results confirm the quality of our predictions: at CaO3?10 C3 , the percentage error compare to theory is round 5%, while at higher Ca is less than 1%.
Finally, we calculate the bubble speed in Newtonian fluid (Fig.6). The bubble moves faster than the liquid flow (U b /U f >1), and its relative velocity tends to increase with the lubricating film thickness. This trend is following the theory (Eq.12), from which our solutions diverge by less than 1%.  The simulated bubble speed in relation with the lubricating film thickness (red dots) compared to the theoretical trend of [19]. In non-Newtonian cases, we change the input parameters of the Carreau model so to investigate the rheological curves shown in Fig.2. Each curve corresponds to a specific Ellis number defined in [23]; in particular, we investigate El=0.1-1-10-100.

Taylor bubble in a shear-thinning fluid
Increasing Ca, the bubble change shape assuming from elliptical to elongated shape (see the well-known from Taylor's experiments [18]); the air tends to expands preferably along the center-axis with a consequent thickening of the lubricating film. When El reduces at fixed Ca, the bubble seems to maintain the equivalent Newtonian shape. However, the lubricating film (Figs.8-9) becomes thinner when the fluid accentuates its shear-thinning behavior. More in general, we observe that the front meniscus assumes a rounded shape accentuating its curvature as Ca increases or El reduces. Such a behavior is coherent to what observed in the simulations with shear-thinning fluids by [26,35].
The region near the tube wall (i.e., where the fluid is at almost-zero velocity) presents a high-viscosity layer that becomes thinner when El decreases (i.e., the fluid is losing the Newtonian behavior, reducing the upper plateau in Fig.2). The thickness of this layer is crucial for the final bubble's shape because it directs the expansion of the air towards the the center-axis of the tube.
Finally, we validate the lubricating film thickness at the lowest Capillary numbers against our recent theory [27]. Results show an excellent agreement with theory, presenting at Ca=10 C3 and 10 C2 an average error by 11% and 3%, respectively.

Conclusion
We presented a numerical analysis of a Taylor  Future developments of this work will assess more in detail the formationof such a high-viscosity layer, investigating its dependence on the rheologyand flow conditions and its e1ects on the lubricating film thickness by theformulation of a scaling law. Further e1orts will investigate the formation ofthe bubble profile in all its parts, as the lubricating