Motion of a confined bubble in a shear-thinning liquid

Abstract The motion of a gaseous Taylor bubble in a capillary tube is typical of many biological and engineering systems, such as small-scale reactors and microfluidic devices. Although the dynamics of a bubble in a Newtonian liquid has been the subject of several studies since the seminal works of Taylor (J. Fluid Mech., vol. 10, issue 2, 1961, pp. 161–165) and Bretherton (J. Fluid Mech., vol. 10, issue 2, 1961, pp. 166–188), the case where the fluid exhibits a shear-thinning behaviour is much less understood. To fill this gap, we study the dynamics of a bubble that moves in a shear-thinning fluid whose viscosity is described by the Ellis viscosity model. With this aim, we derive a lubrication model in the film region to identify the scaling laws for the bubble speed, the film thickness and the pressure drop as a function of the Ellis number and the degree of shear thinning. Our model generalizes Bretherton's results to shear-thinning fluids by identification of a universal scaling law for the effective viscosity that accounts for the interplay of the zero-shear-rate and shear-thinning effects. The film thickness follows a $2/3$ scaling law with respect to the capillary number based on the proposed effective viscosity. The ratio between the bubble speed and the average velocity of the fluid ahead of the bubble is a function of the effective capillary number only. We show that some portions of the bubble are dominated by the zero-shear-rate effect discussing the extent to which the use of the power-law viscosity model can be legitimized. Finally, we study the location of the recirculating vortices ahead of the bubble.

. Sketch of the confined bubble that moves at speed U through a shear-thinning fluid in a channel of half-width R; U ∞ is the average velocity far from the bubble. The shear-thinning fluid and the gaseous bubble are depicted in dark blue and white, respectively. The system of coordinates is placed on the right of the figure just for convenience: in the computation of the rear and front menisci the origin is somewhere in the region CD; AB and EF are the spherical caps, BC and DE are the film regions, while CD is the uniform film thickness region.
A typical set-up that is commonly used to study the dynamics of elongated bubbles is a horizontal capillary tube of radius R, figure 1. In this setting, viscous forces and surface tension dominate over buoyancy and inertia. When the bubble translates at a constant speed, the gas takes a symmetrical bullet shape, commonly called a Taylor bubble, and a thin film of liquid is maintained between the bubble and the tube wall. Knowledge about the thickness of the film and its relationship with the bubble speed are crucial pieces of information for determining the transport rates in practical problems.
In many practical applications the working fluids exhibit a non-Newtonian behaviour (e.g. Savage et al. 2010;Huisman, Friedman & Taborek 2012). Polymers solutions, suspensions, emulsions and biological fluids, just to mention a few, behave like shear-thinning fluids and, therefore, their viscosity is a function of the imposed shear rate. Specifically, at low shear rates, the shear stresses are proportional to the shear rate , where n is the shear-thinning index and Ca * is a modified capillary number for power-law fluids. Unfortunately, the power-law model has an empirical nature and it does not accurately describe the viscosity at small shear rates, leading to large errors when applied to the multiphase scenario (Picchi et al. 2017(Picchi et al. , 2018aPicchi, Ullmann & Brauner 2018b). Both the local velocity field and integral variables (e.g. pressure drop and film thickness) are poorly predicted in free-surface flows due to the unrealistically and unbounded growth of the viscosity at small shear rates. This is confirmed by the work of Hewson, Kapur & Gaskell (2009) that studied the motion of a Taylor bubble through both power-law and an Ellis fluids. The power-law model generates inconsistencies in resolving the low-shear-rate regions of the velocity profile that can be overcome only accounting for the zero-shear-rate effect, see Hewson et al. (2009). Moreira et al. (2020) carried out numerical simulations of Taylor bubbles moving in a Carreau fluid at finite capillary numbers while Kawahara et al. (2015) collected experimental data on a train of Taylor bubbles in shear-thinning polymer solutions. Also, lubricating films of shear-thinning liquids have been the subject of investigation in the context of the drag-out problem, i.e. the Landau & Levich (1942) problem, (Gutfinger & Tallmadge 1965;Tallmadge 1966;Spiers, Subbaraman & Wilkinson 1975;Afanasiev, Münch & Wagner 2007).
However, in all the aforementioned studies, a generalized understanding of the bubble motion that embeds both the low-and high-shear-rate behaviours is still missing, including a generalization of the scaling laws for the film thickness and the bubble speed. To fill this gap, the goal of this paper is to study the dynamics of a confined Taylor bubble that moves in a shear-thinning fluid and to clarify the competition of the zero-shear-rate and the shear-thinning effects on bubble characteristics. We propose a theoretical framework that provides the scaling laws for the film thickness and the bubble speed and quantifies the interplay between low-and high-shear-rate behaviours. Our model generalizes Bretherton's results to shear-thinning fluids by identification of a universal scaling law for the effective viscosity and by proposing a generalization of capillary number that applies to both Newtonian and shear-thinning fluids. We also summarize the limitations of the power-law viscosity model discussing the extent to which its use can be legitimized.
With this aim, we derive a lubrication model for a Taylor bubble that moves in a shear-thinning fluid and whose viscosity is described by the Ellis viscosity model ( § 2). The Ellis model suits our scope since it allows for the derivation of the analytical velocity profile in the approximation of unidirectional free-surface flow. We show that the bubble meniscus obeys a third-order ordinary differential equation for the film thickness ( § 2.3). The model allows for the calculation of the front and rear menisci and the identification of the scaling laws for the film thickness, pressure drop and the bubble speed, see § 3. In § 3.2 we show that, if properly re-scaled, the effective viscosity can be described by a universal law that depends only on the two dimensionless numbers that describes the fluid rheology. Our results converge to Bretherton's theory in the Newtonian limit and we prove that the power-law viscosity model is not an accurate approximation in describing the bubble dynamics. The appearance of recirculating flow patterns ahead and behind the bubble is discussed in § 3.5, while pressure drops across the front and the rear menisci and the interfacial velocity in the film region are presented in Appendices B and C, respectively. The results obtained shed light on the mechanisms that control the motion of a Taylor bubble in a realistic shear-thinning fluid.

Problem formulation
We consider a bubble confined in a horizontal planar channel that advances through a shear-thinning fluid as sketched in figure 1. The bubble is assumed inviscid and translates at a steady velocity U while, far from the bubble, the liquid moves with average velocity U ∞ . We assume that the bubble is sufficiently long so that a region with uniform film thickness h exists and that gravity forces can be neglected. In other words, we assume that the macroscopic Bond number is sufficiently small, Bo = ΔρgR 2 /σ 1 (Δρ and g are the density difference and the gravitational acceleration, respectively). We start from the continuity and Navier-Stokes equations in the shear-thinning fluid without resolving the velocity field in the gas. We look at equations where inertial forces are unimportant in the thin film and the problem is controlled by the competition between viscous and surface tension forces only. The governing equations in the liquid are where p and τ are the pressure and the shear stress tensor; u = (u, v) is the velocity vector with u and v representing the velocity components in the x and y directions, respectively. The boundary conditions are no slip and no penetration at the channel wall and a free shear-stress condition at the bubble-fluid interface.
The shear-thinning fluid is assumed to be an incompressible generalized Newtonian fluid of density ρ and effective viscosity μ, whose constitutive relation can be formulated as whereγ = ∇u + (∇u) T is the rate-of-strain tensor. Since we are interested in understanding the interplay between the shear-thinning and the zero-shear-rate effects on the bubble dynamics, we chose the Ellis viscosity model (Matsuhisa & Bird 1965;Reiner 1965;Bird et al. 1987). The constitutive law of an Ellis fluid is given by where τ = √ 1/2(τ : τ ) is the effective shear stress (or the magnitude of the shear-stress tensor). The model converges to Newtonian viscosity μ 0 , the zero-shear-rate viscosity, at small shear rates. The index α represents the degree of shear thinning and it assumes values α ∈ [1, ∞); the greater the value of α, the greater is the extent of shear thinning. The constant τ 1/2 controls the onset of the shear-thinning effect representing the effective shear stress at which the viscosity is 50 % of the Newtonian limit. The greater the value of τ 1/2 , the greater is the range of shear rates with the Newtonian viscosity. The Ellis model describes well the viscosity of many aqueous solutions (e.g. carboxymethyl cellulose and Carbopol solutions), where the typical values of the rheological constants are in the range α ∈ (1, 3) and τ 1/2 = O(10 −3 ÷ 10), and the viscosity of many polymers, where α ∈ (1, 4) and τ 1/2 = O(1 ÷ 10), see Bird et al. (1987).
In the next section, we derive a lubrication model for a Taylor bubble when it reaches a steady state and moves along the tube with constant speed U. Specifically, once formed, the uniform film thickness h is constant with time in a reference frame that moves at the bubble speed (x − Ut). This assumption is supported by the numerical simulations of Moreira et al. (2020) and the experiments of Kawahara et al. (2015). Alternatively, in problems where the main interest is in characterizing the transient evolution of the bubble, one should derive an evolution equation for the film thickness in time (see for example Zhang & Lister 1999;Eggers & Fontelos 2015;Garg et al. 2017). In addition, we assume that the film thickness remains high enough to justify neglecting the effect of long-range intermolecular forces (e.g. van der Waals, disjoining pressure). Therefore, the possible rupture of the film and dewetting of the channel wall due to those destabilizing effects (see Zhang & Lister 1999;Diez & Kondic 2007;Garg et al. 2017) are out of the scope of this work.

Lubrication model
We consider the portions BC and DE of the bubble in figure 1, where the curvature of the meniscus is established and the interface slope is small, |dy i /dx| 1, with y i denoting the location of the interface. The idea is to simplify the governing equations by introducing a small scaling parameter ε, defined as the ratio of the uniform film thickness h to the unknown length of the film region along the flow direction, h/ε. This approach is typical of thin-film lubrication models (Eggers & Fontelos 2015) and, in the limit of ε 1, it allows for the identification of the scaling laws in Landau-Levich-Derjaguin-Bretherton problems (see de Gennes, Brochard-Wyart & Quéré 2003; Stone 2010).
Assuming that the velocity in the axial direction scales with the bubble speed U, the following dimensionless variables are introduced where V is the characteristic scale of velocity in the y direction. The continuity equation (2.1a) suggests that V = εU, and we chose the zero-shear-rate viscosity μ 0 as representative scale for the effective viscosity. The pressure in the film region scales with the normal pressure that can be estimated as p ≈ σ (d 2 y i /dx 2 ) ∼ σ ε 2 /h. Normalizing the momentum equations (2.1b) using (2.5) yields are the Reynolds and the capillary numbers both based on the zero-shear-rate viscosity.
Assuming that the scaling parameter is small, ε 1, and that inertial terms can be neglected, εRe 1, (2.6) simplifies to (2.9) Equation (2.9) shows that the viscous stress balances the axial pressure gradient only if ε 3 ≈ Ca and, therefore, thex variable can be expressed in terms of the capillary number asx = x/hCa −1/3 . This allows simplification of (2.7) to ∂p/∂ỹ = 0, whereby the pressure is a function of the axial variable only. Note that the model is not applicable near the cap of the bubble where the curvature scales with the channel half-width R, and the slope of the interface is not small enough to justify the lubrication approximation. The effective shear stress reduces to which is equal to the dimensionless shear stress in the film,τ xy ≈μ∂ũ/∂ỹ. The corresponding effective viscosity can be written as where El is the Ellis number. The Ellis number can be interpreted as the ratio between the characteristics shear rate of the fluid, τ 1/2 /μ 0 , and the representative shear rate in the film, U/h. When El → ∞ the onset of the shear-thinning effect is delayed to an infinite shear rate and the model reduces to a Newtonian fluid of viscosity μ 0 (andμ = 1), figure 2(a).
We refer to this case as the Newtonian limit. Instead, when the shear rate is sufficiently high, or El is small such that the Newtonian plateau shrinks to extremely low shear rates, the shear-thinning effects dominates in (2.11) and the viscosity matches that of a power-law 10 -2 Power-law Power-law 10 -4 10 -6 10 -3 10 -2 10 -1 10 0 10 0 10 2 10 -2 10 -4 10 -6 10 0 μ 10 -3 10 -2 10 -1 10 0 fluid with a slope (1 − α)/α, figure 2. As α → 1 the fluid exhibits a power-law behaviour almost in the entire range of shear rates, while for α = 1,μ = 1/2. At the interface we ensure the continuity of normal and tangential stresses and the boundary conditions read Since Ca ≈ ε 3 and ε 1, (2.12) further simplifies to (2.14) The assumption that ε 1 implies that the slope of the interface in the film region is small, namely dη/dx Ca −1/3 . A formulation of the boundary condition that relaxes this constraint is discussed in Park & Homsy (1984) and Reinelt (1987).

Film equation
We integrate (2.9) subjected to (2.13) atỹ = η and the no-slip condition atỹ = 0 to obtain the (instantaneous) velocity profile of shear-thinning fluid The velocity profile is composed by a Newtonian and a shear-thinning parts, reflecting the structure of the Ellis rheological model that can be seen as the sum of a Newtonian and a power-law viscosity model, see Gutfinger & Tallmadge (1965). The details concerning the derivation of the velocity profile are listed in Appendix A. Then, we compute the volume flux by integrating the velocity profile (2.15) over the film thickness Looking at the problem from a reference frame moving with the bubble at the speed U, the mass balance between the region of uniform thickness CD and the film region CB (see In the region CD the liquid is at rest since η = 1 and the pressure is constant from (2.14).
In terms of dimensionless coordinates, (2.17) becomes Combining (2.16) and (2.18) and using the boundary condition (2.14) to express the driving force, dp/dx = −d 3 η/dx 3 , we obtain we finally get the following ordinary differential equation for η where the bubble profile η is a function of ξ , α and El only. Equation (2.20) describes the meniscus between the uniform film and the caps at both the bubble ends in a moving reference frame (x − Ut). The left-hand side of (2.20) is composed of a Newtonian and a shear-thinning term, indicated as I and II, respectively. When the II becomes negligible, (2.20) converges to the classical result obtained by Bretherton (1961) for a Newtonian liquid, We recover the Newtonian liquid either for El → ∞, or for specific combinations of El and α so that the shear-thinning effect vanishes, see figure 2. Also, when the third derivative goes to zero, d 3 η/dξ 3 → 0, (2.20) reduces to (2.21), as it will be discussed in detail in § 3.1. If we artificially remove the term I from (2.20), we recover the equation for a power-law fluid  Kamişli & Ryan (1999, 2001 for power-law fluids. Although mathematically correct, the model for power-law fluid leads to non-physical results when integrated starting from the uniform thickness region and the extent to which its use could be legitimized will be discussed in § 3.1. The general solution of (2.20) can be found numerically starting from the region CD and integrating in the directions of the front and rear menisci. To do that, it is convenient to discuss few limiting cases where the solution can be obtained analytically ( § 2.4).

Scaling of the film equation
In regions near C and near D the film thickness approaches a uniform value (y i h) and we can write η ≈ 1 + η 1 + · · · with η 1 1. (2.23) Upon substituting (2.23) into (2.20) we obtain ( 2.24) Expanding the two binomials in (2.24) into Taylor series, truncating after the second term and recalling that η 1 = η − 1, we obtain an approximation of the film equation near the uniform film and it admits an analytical solution where a, b, c are constants. This is the same approximation valid for Newtonian fluids and it holds also for an Ellis fluid since the third derivative in this region is small, d 3 η/dξ 3 1; this assumption will be justified a posteriori in § 3.1. In the following, we refer to portions of the bubble where (2.27) applies as the exponential regions, see figure 3.
In regions where η 1 and the slope of the interface is small enough, dy i /dx 1 (or, in terms of dimensionless variables dη/dξ Ca −1/3 ), so that the lubrication approximation holds, (2.20) reduces to d 3 η dξ 3 ≈ 0 with η 1. (2.28) Integration yields a parabolic behaviour with respect to ξ where P, W and Z are constants to be determined from the numerical solution of (2.20), see § 2.5. Specifically, the coefficient P is the dimensionless curvature of the meniscus when η 1, i.e. the second derivative of η with respect to ξ , P = d 2 η/dξ 2 . Given El and α, the value of P can be obtained from the numerical solution far from the uniform film region for both the front and the rear menisci. The determination of the coefficients W and Z requires numerical fitting: W is made zero by shifting the solution along ξ , while Z is fitted from a window that gives, in the Newtonian limit, the same value obtained by Bretherton (1961) (Z = 2.79 in the front and Z = −0.72 in the rear). The portions of the bubble where (2.29) applies are referred as the parabolic regions, see figure 3.

Boundary conditions and numerical solution
The film shape between the uniform film and the caps at both ends is obtained by solving (2.20) using the fully implicit solver ode15i of Matlab. The front and the rear menisci are treated separately integrating (2.20) with a different set of boundary conditions. The front meniscus is obtained by integrating (2.20) starting from somewhere near C (figure 1) towards positive ξ . We chose to set the origin at ξ = 0, assuming that the uniform film region extends to ξ → −∞, or, alternatively, that η(−∞) = 1. This means that, near C, we can approximate the solution with the exponential approximation (2.27). In the front, the cosine and sine terms in (2.27) are small due to the damping effect of the negative exponentials, and (2.27) can be simplified as (2.30) For specified values of El and α, the solution is unique. Starting the integration from ξ = 0 we set a = 10 −5 and Different values of a result only in a shift in the origin of ξ . The (translated) front profile is shown in figure 3(a). The exponential solution is also plotted indicating that, as ξ increases, the numerical solution diverges from (2.30).
The profile in the rear meniscus is obtained by integrating (2.20) towards negative ξ starting from the boundary condition η(∞) = 1. Near D, we can approximate the solution with the exponential function (2.27) that, for negative ξ , can be simplified neglecting the positive exponential as (2.32) Choosing ξ = 0 as the origin we get In (2.33a-d) there are two disposable constants: b has the effect of shifting the origin of the bubble profile only and we arbitrarily set the value to b = 10 −4 . The other constant determines the curvature of the menisci at η 1, i.e. for ξ → −∞. It is a common choice to assume that the curvature in the rear matches the one in the front, see Bretherton (1961), imposing (2.34) The boundary condition is satisfied by an iterative loop on the constant c in (2.33a-d); this loop imposes that the second derivative -the curvature -of η in the parabolic region equals the one of the front meniscus calculated for given values of El and α. A plot of the calculated rear meniscus is given in 3(b). Differently from the front meniscus that is monotonic, the rear meniscus develops the typical oscillations. A discussion about the relaxation of the boundary condition (2.34) is given in Appendix D.

Results and discussion
3.1. Characteristics of the front meniscus We compute the front meniscus as described in § 2.5 aiming at exploring the effect of the rheology of the shear-thinning liquid on bubble characteristics. Differently from the Newtonian case, where the profile is uniquely determined by the dimensionless variables η and ξ , here, the profile also depends on El and α.
The Ellis number controls the onset of the shear-thinning effect. When El → ∞, the front meniscus converges to the Newtonian limit and the problem does not depend anymore on α, figure 4(a). Instead, as El decreases, the Newtonian plateau in the viscosity curve shrinks, see figure 2(a), and the dimensionless curvature of the meniscus d 2 η/dξ 2 when η 1 decreases, see figure 4(b). The change in the meniscus shape is an effect of the local variation in the effective viscosity of the liquid, as demonstrated in figure 5. The effective viscosityμ| w is computed using (2.11) based on the shear rate at the channel wall (ỹ = 0) defined as dũ dỹ w = − dp dx η − El 1−α dp dx dp dx α−1 η α . (3.1) In the vicinity of the uniform thickness region, the effective viscosity isμ| w = 1 since the fluid is at rest. This means that we can always identify a region of the bubble where the liquid exhibits a Newtonian behaviour. The existence of a high (Newtonian) viscosity  region within the uniform film thickness is indeed confirmed by the viscosity field obtained from numerical simulations by Moreira et al. (2020). To quantify the shear-thinning effect, the weights of the Newtonian and shear-thinning terms are plotted in figure 5 as where I and II refers to (2.20). For El = 10 the shear-thinning term is negligible over the entire meniscus and the profile is almost indistinguishable from a Newtonian profile. Instead, for El = 1, 10 −1 , 10 −5 we can identify a region of the bubble where the shear-thinning term in (2.20) is comparable to the Newtonian one, figure 5(c), or entirely dominates the solution, figure 5(d). We define the shear-thinning region (depicted in grey) when I N ≤ 0.95. The terms (3.2a,b) suggest that the shear-thinning effect is important only when the film starts growing rapidly before entering the parabolic region. There, the shear rate has a maximum that corresponds to a minimum in the effective viscosity. As expected, the lower El is, the lower is the minimum ofμ| w . Far from the uniform thickness region where η is large, the third derivative d 3 η/dξ 3 → 0 and, as a consequence, the effective viscosity reduces rapidly, see (3.1). Although this would imply the tendency of reaching the Newtonian limit again, we remind that the model is not applicable too far from the region CD.
The fact that the shear-thinning effect arises only in a small portion of the bubble indicates that the use of the power-law model for describing the meniscus is physically incorrect. In particular, since in the vicinity of the uniform film region, the problem is dominated by the Newtonian term (figure 5), the integration of (2.22) would lead to non-physical results. The power-law effective viscosity would grow to the infinite near the origin since dũ/dỹ ≈ 0. The choice made by Kamişli & Ryan (1999) of solving (2.22) starting from a Newtonian initial condition, but not accounting for the Newtonian term, does not seem a correct representation of the physics of the problem. The shear-thinning effect also delays the transition to the parabolic region that is characterized by a constant dimensionless curvature and d 3 η/dξ 3 ≈ 0, see § 2.4. As El decreases, the starting point of the parabolic region shifts to higher ξ , as can be seen by analysis of the rate at which the third derivative approaches a zero value. In the Newtonian limit, the third derivative in (2.21) decays proportionally to 1/η 2 . In case the shear-thinning effect dominates, such as for El = 10 −5 in figure 5(d), the third derivative in (2.22) decays proportionally to 1/η (1+α)/α , (3.4) Since α > 1, the onset of the parabolic region is delayed. Changes in the degree of shear thinning α result in two different behaviours. When the shear rate is smaller than 0.2 (the value at whichμ = 1/2), the system tends to the Newtonian limit for α → ∞. Instead, when the shear rate is greater than 0.2, an increase in α translates in a decrease of the effective viscosity, see figure 2(b). Accordingly, in figure 6, the front meniscus approaches the Newtonian limit for El = 1 while it shows the opposite behaviour for El = 10 −2 . Since the effective shear rate of the problem scales as γ e ∼ 1/El, the two opposing behaviours depend on whether the Ellis number is greater or smaller of a critical value. For α = 1 we recover a Newtonian fluid withμ = 1/2 and the meniscus is identical to the one obtained in the Newtonian limit with a change of variables ξ * = ξ/ 3 √ 2. The characteristics of the front meniscus can be summarized looking at the scaling behaviour of the coefficients in the parabolic solution (2.29). Figure 7(a) shows P as a function of El and α. When El → ∞ all the curves approach the Newtonian limit and P = 0.643 as obtained by Bretherton (1961). The curvature of the Newtonian limit is the maximal curvature that the meniscus can assume. As El → 0, instead, the curvature decreases significantly and, when α → 1, the shear-thinning effect delays the transition to the Newtonian limit to much higher El. Figure 7(b) shows the trends for the coefficient Z in (2.29).

Scaling of the film thickness and generalization of the capillary number
The model for the film thickness is obtained by matching the curvature of the parabolic region near B with the curvature of the spherical cap AB, see figure 1. Specifically, near B the curvature is constant d 2 y i  . Coefficients P and Z in (2.29) that characterize the front meniscus in the parabolic region as a function of the Ellis number and the degree of shear thinning α. Note that the coefficient P is also the dimensionless curvature of the meniscus for η 1, P = d 2 η/dξ 2 .
the Newtonian limit. The scaling law for h/R is summarized in figure 8. The film thickness decreases with decreasing of the Ellis number since the effective viscosity of the problem diminishes, figure 8(a). Therefore, when the shear-thinning effect dominates (i.e. low El), the bubble forms a thinner film compared with a Newtonian fluid with the same capillary number. Variations in the degree of shear thinning α produce two opposite behaviours depending on whether the effective shear rate (or 1/El) exceeds the critical value, see 8(b,c). The model for the film thickness (3.6) is still based on the definition of the capillary number defined with the zero-shear-rate viscosity, see (2.8a,b). It would be desirable to obtain a formulation of the capillary number that embeds both the zero-shear-rate and the shear-thinning effects. With this aim, we can recast (3.6) in the following form introducing the effective capillary number Ca e h R = 0.643(3 Ca e ) 2/3 with Ca e = μ e U σ , where μ e is the effective viscosity defined as (3.8) A plot of the effective viscosity as a function of the effective shear rate 1/El is provided in figure 9(a). Since the determination of P requires the numerical solution of (2.20), we use the numerical results to get a fitting curve for the effective viscosity. The data of figure 9(a) collapse over the same master curve, figure 9(b), if we define  et al. (2020). In particular, the film thickness decreases with increasing of the degree of shear thinning at a fixed capillary number. Unfortunately, their data cannot be used to validate our model since the flow rates necessary to define the capillary number are not reported in the manuscript. Also, the data collected by Kamişli & Ryan (1999) and by Kawahara et al. (2015) cannot be considered to validate the scaling law for the film thickness because the rheological data do not include the low-shear-rate behaviour and only the power-law behaviour has been provided. For a complete validation of our model, more experiments or simulations performed with realistic shear-thinning fluids would be valuable and are encouraged.
3.3. Bubble speed A measure of the ratio between the bubble speed U and the average velocity of the fluid ahead of the bubble U ∞ is a critical parameter in the design of many applications that involve bubbles and trains of bubbles. We can derive the scaling law for U/U ∞ applying the mass balance to the problem. Specifically, if we write the balance between the open tube and the uniform film region in a reference frame that moves with the bubble speed U (3.12) Combining (3.12) with (3.6) or (3.7) we obtain the scaling law valid for the two-plate geometry The bubble always flows faster compared with the fluid ahead from the bubble, U/U ∞ > 1, and, only in the limit of Ca → 0, the speed of the bubble is almost indistinguishable from U ∞ , see figure 10. In general, when the shear-thinning effect plays a role (i.e. low El), the bubble flows slowly compared with a Newtonian bubble with the same capillary number Ca. The lower the effective viscosity is, the slower is the bubble compared with U ∞ . The effect of α is summarized in figure 10(b,c), where we see that increasing the degree of shear thinning can increase the bubble velocity at El = 1. As expected the scaling law for U/U ∞ converges to the one obtained by Bretherton (1961) for El → ∞. The trends of figure 10 are in a qualitative agreement with the simulations of Moreira et al. (2020) where the bubble speed decreases with increasing of the degree of shear thinning at a fixed capillary number. Using the definition of the generalized capillary number Ca e defined in § 3.2, the ratio U/U ∞ is monotonic with Ca e and all the curves collapse on the same one depicted in 10(d). This analysis can be generalized to the case of a pipe geometry by referring to the film holdup instead of the film thickness in (3.11). In the case of circular geometry, the film holdup scales with (3.14) In figure 10(d) we can see that the speed ratio U/U ∞ attains a higher value in a pipe compared with the two-plate geometry.

Characteristics of rear meniscus and limitations of the lubrication approximation
We compute the profile of the rear meniscus as described in § 2.5. In terms of the curvature far from the uniform film region, the rear meniscus resembles the same trends of the bubble front, see figure 11. The main difference is the presence of oscillations in the bubble profile that stretches along ξ as the shear-thinning effect becomes more important. As the curvature P decreases, the minimum of the profile shifts closer to the channel wall. I N and I PL show that, similarly to the front, the region in the vicinity of the uniform film thickness is dominated by the Newtonian term. Instead, there is a competition between the two near the minimum of the profile, while at ξ → −∞, the shear-thinning term decays and the problem converges again to the Newtonian one. Interestingly, the effective viscosity in the rear is not a monotonic function of ξ . Due to oscillations in the profile, the pressure gradient (dp/dx = −d 3 η/dξ 3 ) changes sign and thereby inverting locally the direction of the flux. The effective viscosity reflects this behaviour showing a certain number of spikes in correspondence of such changes of sign. This is an apparent inconsistency of the model that can be attributed to the lubrication approximation. In fact, variations in the ξ direction are accounted only as a first-order approximation in the scale parameter ε, but, in the vicinity of the oscillations, the assumption that ε 1 should be relaxed. The effective viscosityμ w presented in figure 12 is the effective viscosity calculated based on the shear rate at the wall neglecting the effect in the ξ -directionγ (3.15) Equation (3.15) converges to (3.1) when ε 1. However, in the vicinity of the oscillations the term ∂ũ/∂ξ is not of the order one and the effective shear rate needs to be corrected. A way for estimating the velocity gradient in the axial direction could be to consider the gradient of the average velocity in the film ∂ũ/∂ξ ≈ ∂Ũ av /∂ξ yielding where β is a parameter that gives the order of magnitude of the axial derivative term. In our analysis, we do not have a way for estimating precisely the magnitude of the axial derivative, but in figure 13 we show the viscosity profile regularized for different values of β. Even a small β is sufficient to smooth the viscosity profile. To conclude, differently from the front meniscus where the viscosity profile is regular, a correction that accounts for the axial derivative of the velocity in the computation of the effective shear rate is necessary for a correct interpretation of the rear viscosity field.
3.5. Recirculating flow patterns In Taylor bubble flow, recirculating flow regions can form ahead of and behind the bubble. In the low capillary number limit, in fact, the bubble is surrounded by a thin film and a qualitative sketch of the streamlines in a reference frame attached to the bubble is provided in figure 14. As long as the film thickness is small, h/R 1, and the bubble flows slower compared with the maximum liquid velocity ahead, an external re-circulation flow pattern (centred at the location y 0 ) is present. However, when the film thickness is greater than a critical value, the bubble can flow faster than maximum liquid velocity ahead and, therefore, the external re-circulation regions are not present, as discussed by Taylor (1961), Giavedoni & Saita (1997), Thulasidas, Abraham & Cerro (1997, Rocha, Miranda & Campos (2017) and Balestra et al. (2018) for Newtonian fluids. Since a fundamental understanding of the location of the vortices is a crucial piece of information when trains of Taylor bubbles in capillaries are used to enhance heat and mass transfer, we compute the critical film thickness for the appearance of flow re-circulation and the location of the vortices for the case of shear-thinning fluids. In a coordinate system attached to the bubble, a recirculating pattern is expected when the maximal velocity ahead u max ∞ exceeds the bubble velocity U, i.e. u max ∞ /U > 1 in figure 1. Assuming a fully developed flow ahead of the bubble with average velocity U ∞ , the velocity profile reads are dimensionless variables based on the channel half-width R and U ∞ . The pressure gradient in (3.17) is determined numerically satisfying the flow rate constraint The maximum of the liquid velocity is at the channel centre (i.e. Y = 1) (3.20) Differently from the case of Newtonian fluids, where the ratio between the maximal and the average velocity is 3/2, the shear-thinning effect reduces the ratio u max ∞ /U ∞ due to the 918 A7-21  flattening of the velocity profile in the vicinity of the channel axis, see figure 15(a). The condition for the appearance of recirculating patterns is u max and, using (3.12) in (3.21), we can identify the critical film thickness as h * R <ũ max − 1 u max . (3.22) The critical thickness for the appearance of the recirculating patterns depends on El ∞ and α as shown in figure 15(b). In the Newtonian limit (for El ∞ → ∞), h * /R = 1/3 in agreement with Balestra et al. (2018). The shear-thinning effect reduces the critical film thickness at small El ∞ , in particular when α is large. The centre of the recirculation zone y 0 can be estimated by looking for the point where the velocity equals the bubble velocity. In terms of dimensionless coordinates, Y 0 = y 0 /R is found by satisfyingũ ∞ (Y 0 ) = U/U ∞ combined with (3.12). Figure 16(a) shows that, as the shear-thinning effect becomes important, Y 0 shifts closer to the channel wall (for h/R 1). We can also estimate the location of the dividing streamline separating the circulating vortex and the liquid flowing towards the film. The location of the dividing streamline, y d in figure 14, is computed knowing that the liquid flow rate from y = 0 to y = y d equals the flow rate in the uniform film region, see Thulasidas et al. (1997) and Picchi et al. (2018b). In fact, the net flow rate calculated over the vortex region (depicted in red in figure 14) is zero in the moving reference frame. The mass balance between the uniform film and the portion of the fluid with non-zero flow rate yields (3.23) Casting (3.23) in terms of dimensionless variables and using (3.12) we obtain an equation El ∞ = 10 -3 El ∞ = 10 -2 El ∞ = 10 -1 El ∞ = 1 El ∞ = 10 Newtonian Solving (3.24), we obtain that the location of the dividing streamlines scales with the dimensionless film thickness, Y d ∼ h/R, and it is practically independent on the shear-thinning effect as long as h/R 1, figure 16(b). The fact that the recirculating vortices extend close to the channel wall in the simulations of both Newtonian and shear-thinning fluids (Balestra et al. 2018;Moreira et al. 2020) suggests that such scaling could be an universal behaviour. The same arguments proposed in this section also apply to the region far behind the bubble.

Conclusions
In this paper, we study the motion of a confined Taylor bubble in a shear-thinning fluid. We derive an equation that describes the meniscus profile for an Ellis fluid accounting for both zero-shear-rate and shear-thinning effects. Our analysis identifies a universal scaling law for the generalized effective viscosity of the problem that applies to both Newtonian and shear-thinning fluids, and allows identifying the scaling laws for the bubble velocity, the film thickness and the pressure drop. In particular, the two-third scaling law for the film thickness holds only if the capillary number is formulated based on the generalized effective viscosity, h/R ∼ Ca 2 e /3. We show that some portions of the bubble are dominated by zero-shear-rate effects (e.g. the vicinity of the uniform film thickness), while near the nose of the bubble, it is difficult to discern between the two effects. To overcome such difficulties, we quantify clearly the relative importance of the shear-thinning effect (and the consequent reduction in the effective viscosity) over the meniscus. Our study also clarifies the limitations of the power-law viscosity model showing that its use in the Taylor bubble problem leads to non-physical results. A full rheological characterization, including the low-shear-rate behaviour of the working fluids, is encouraged in future experiments so that new data can be used for future validation and extension of the model.
Despite the fact that the motivation of our work is oriented to microfluidic applications, the scaling relations obtained for the film thickness may serve as a generalization of the well-known Landau-Levich-Derjaguin-Bretherton problem to the case of inelastic shear-thinning fluids. As shown by Stone (2010), there are several practical problems that involve thin films with a lot of similarities from a fluid mechanics perspectives, such as many classes of coating problems (e.g. coating a plate by vertical withdrawal from a bath of liquid or roll coating of a horizontal rotating cylinder partially immersed in a bath of liquid). Since in all these applications the film thickness is not known a priori, the identification of such scaling relations is of great importance. Furthermore, the study of the speed of an isolated Taylor bubble through a shear-thinning fluid can be used to construct more sophisticated models for trains of bubbles (e.g. slug flows).
Declaration of interests. The authors report no conflict of interest.
D. Picchi https://orcid.org/0000-0002-2619-104X. region and stabilizes at η 1, figure 17(a). The limiting value is For a correct interpretation of figure 17, recall that the theoretical limit of 1/2 has a physical meaning only if the film thickness is small compared with the pipe radius, ηh/R 1.

Appendix C. Pressure drop
The pressure drop across the bubble can be estimated following the same approach by Bretherton (1961). The idea is to approximate the bubble profile in the film region with an effective profile that has constant curvature, so that the Young-Laplace law can be used. In the front, upon substituting (3.6) into (2.29) we obtain a continuation of the spherical region through the region BC, figure 1, Equation (C1) has a tangent parallel to the wall at x = 0. We use the distance to the wall at x = 0 to compute the approximated curvature with Ca 1 (a Taylor expansion is used to simplify the formula of the curvature). Using the Young-Laplace law, the dimensionless pressure drop in the front meniscus yields where Δp originates both from the surface tension stresses across the interface (static pressure drop) and viscous stresses due to the flow in the film region (dynamic where Z r is the fitting parameter obtained by shifting the solution in order to get W = 0 in (2.29) and Z r = −0.72 in the Newtonian limit (the same as in Bretherton 1961). The pressure drop in the rear has an opposite sign with respect to the one at the front, see figure 19(a) for the case with α = 2 and figure 19(b) for the case with El = 10 −1 .