Flows inside polymer microfluidic droplets: Role of elasticity

,


Introduction
Microdroplets are widely employed in chemical, material, biological or soft matter applications, due to their numerous advantages such as small sample volumes, accurate manipulation high throughput, integration miniaturization and two-phase interface cross-contamination avoidance (Seemann et al., 2012;Chen et al., 2022;Zhou et al., 2020;Li et al., 2021;Zeng and Fu, 2020).In many microdroplet applications (Ma et al., 2021;Marshall and Walker, 2019;Samandari et al., 2021) the fluids can be viscoelastic, due to the presence of polymers or macromolecules, which may lead to complex hydrodynamic phenomena.Viscoelasticity affects the formation of droplets (Vagner et al., 2021;Yu et al., 2020;Xue et al., 2020), droplet dynamics (Asghari et al., 2020;Hatch et al., 2013;Smith and Bertola, 2010), as well as the flow fields inside the droplets.Elucidating the impact of viscoelasticity on the flow behavior inside droplets can not only improve our fundamental understanding of non-Newtonian droplet flow dynamics in microchannels (Asghari et al., 2020;Hatch et al., 2013), but also provide tools to control the mixing or efficiency of droplet reactions involving complex fluids and the mechanical microenvironment around encapsulated cells (Chen et al., 2018;Ma et al., 2015).
The flow patterns inside Newtonian droplets have been the focus of extensive research.They are typically characterised by the existence of two counter rotating vortices.However, a number of flow transitions have been observed as a result of the competitive action of two-phase viscous shear stresses, wall confinement and interfacial tension, giving rise to more complex flow topologies.The capillary number Ca-defined as Ca = μ c V c /γ, where μ c , V c and γ are the continuous phase viscosity, velocity and interfacial tension respectively-is one of the key parameters determining Newtonian microdroplet flows.In rectangular microchannels-unlike cylindrical capillaries-the droplet flow topology is also influenced by the existence of gutter flows at the corners of the channel (Jakiela et al., 2012;Li et al., 2021;Li et al., 2020), which also results in droplet velocities to deviate significantly from the average flow velocity-by up to 30% (Helmers et al., 2019;Jose and Cubaud, 2014).By means of micro particle image velocimetry (micro-PIV), Kinoshita et al. (Kinoshita et al., 2007) resolved the three dimensional flow structure inside droplets in rectangular channels; they observed a double eddy pattern which they attributed to the role of the continuous phase shear.
Jakiela et al. (Jakiela et al., 2012) experimentally observed a discontinuous flow transition from a double to a six-eddy pattern inside microdroplets in rectangular channels, at Ca ~ 10 -3 , due to additional smaller eddies forming at the cups of the droplet.This transition is opposite to one theoretically observed by Hodges (Hodges et al., 2004) for viscous droplets in capillaries.Similar transitions to Hodges have been reported by Liu et al. (Liu et al., 2017) using sunflower oil as the continuous phase, which is more viscous than the hexadecane used by Jakiela (Jakiela et al., 2012), resulting thus in higher viscous stresses in the continuous phase.Besides Ca, Ma et al. (Ma et al., 2014) also studied the effects of viscosity ratio, droplet size, and interfacial tension on the flow fields inside droplets, and found that the flow field depends heavily on the viscosity contrast between the two-phases rather than the Marangoni effect, whereas the Ca and droplet size were found to have no effect on the flow topology.
Kovalev et al. (Kovalev et al., 2018) studied the flow fields inside droplets of various shapes moving in a highly viscous oil and extremely low viscosity ratios; they observed a variety of flow topologies ranging from the typical structure of two counter rotating vortices to one in which no vortices were present.These were attributed to the different plug shapes and described by means of a modified capillary number multiplied by the ratio of dispersed to continuous phase flowrates.Kovalchuk et al. (Kovalchuk and Simmons, 2021) explored the effect of surfactants and found that the Marangoni stresses accelerated the gutter flow.Rao and Wong (Rao and Wong, 2018) found that the continuous phase exhibits different motions in the gutter flow that influence internal flows.Our previous work indicated that the competition between oil film resistance and corner flow determines whether the six-eddy structure appears inside droplets (Li et al., 2020).Mießner et al. (Mießner et al., 2020) simultaneously measured the flow field inside and outside droplets moving in rectangular microchannels; they reconstructed the 3D dimensional flow field from 2D measurements in various planes and were able to describe the flow features in great detail; the strength of main and secondary vortices in the rear and front caps as well as the vortices outside of the droplet were quantified and the stagnation points in the flow investigated.Helmers et al. (Helmers et al., 2022) extended that work in an attempt to characterise topological transitions in droplet flows as a function of Reynolds number and viscosity contract λ.The location of the main vortex core and its diameter varied with increasing the combined parameter Re, λ and the changes were attributed to the influence of two momentum transfer mechanisms; a wall film based one and transfer across the interface in the gutter flows.
The flows inside non-Newtonian droplets have received relatively less attention.The majority of the published literature focusses on the effects of shear thinning rheology on droplet formation/characteristics, the flow inside but also outside droplets.Roumpea et al. (Roumpea et al., 2017) studied experimentally the flows in oil plugs dispersed in a shear thinning fluid using two colour PIV; the shear thinning affected the shape of the profiles in the continuous phase as well as the circulation times.Shear thinning has been found to affect the formation of microdroplets due to the shear varying viscosity (Gu and Liow, 2011;Wong et al., 2017;Vagner et al., 2018).In our previous work (Li et al., 2021) we found that the flow structure inside shear thinning droplets was highly depended on the extent of shear thinning of the fluid; we also observed teardrop-shaped droplets forming at high flow rates due for strongly shear thinning droplets.
In most of these non-Newtonian studies the viscoelasticity of the working fluids is not well characterized and it is often difficult to decouple the effects of elasticity and shear thinning.The latter can be achieved using Boger fluids which have a constant shear viscosity but are elastic.Fluid elasticity can induce secondary flows and instabilities in micro and macroscale flows, both in the absence (purely elastic) (Poole, 2019;Qin et al., 2020) or presence of inertia (elastoinertia instabilities) (Burshtein et al., 2017;Choueiri et al., 2021;Lacassagne et al., 2020;Lacassagne et al., 2021;Boulafentis et al., 2023).
Boger fluids have been employed in interfacial microfluidic studies; however, the focus has been primarily on filament thinning and droplet formation (Vagner et al., 2021;Yu et al., 2020;Xue et al., 2020) where elasticity leads to the formation of threads and beads.Recent studies with viscoelastic continuous phases have shown that elastic instabilities can induce droplet oscillations and their permanent trapping in converging diverging microchannels or mobilize capillary entrapments (Xie et al., 2022;Shakeri et al., 2021).Experiments with macroscale viscoelastic droplets falling into a viscous fluid have shown fascinating topologies in which the elastic stresses induce the formation of a filament of the continuous phase inside the droplet and eventually toroidal droplets (Sostarecz and Belmonte, 2003).However, the role of elasticity on the flow behavior and topology inside microdroplets in rectangular channels has not been reported yet, to the best of our knowledge.
The lack of systematic characterisation of viscoelastic microdroplets and the wealth of elasticity driven phenomena have thus motivated the present study.Using Boger fluids and micro-PIV we study the flow characteristics inside viscoelastic droplets for a range of Weissenberg (Wi), capillary (Ca) and elasticity (El) numbers.The results are compared against those for inelastic Newtonian droplets and discussed in the light of the relative importance of viscous and elastic forces.

Materials
Table 1 summarizes the fluids employed in the present study and their physical properties.Sunflower oil was used as the continuous phase due to its high viscosity compared to the aqueous phase and no swelling effects on Polydimethylsiloxane (PDMS, Sylgard 184, Dow Corning Co.) (Jakiela et al., 2012;Li et al., 2021;Li et al., 2020;Liu et al., 2017).0.15 wt% Span 80 (Sigma-Aldrich, UK) was added to the oil phase to reduce interfacial tension and promote droplet formation.Water-glycerol (G9012, Sigma-Aldrich, UK) solutions were used as the dispersed phase.Polyacrylamide (PAAM, molecular weight M w = 5.5 × 10 6 Da, Sigma-Aldrich, UK) was added at two different concentrations, 250 ppm and 500 ppm, to create the Boger droplets.The critical concentration C* of the polymer solutions, estimated from the Solomon method (Solomon and Ciutǎ, 1962), was found to be around 530 ~ 560 ppm, hence, we can assume that the dispersed polymer phases employed in the present work are in the dilute regime.The extensional rheological properties of the PAAM solutions were characterized though capillary break up experiments as shown in Fig. 1.The elasticity of the fluid becomes evident; in the absence of polymers, the filament breaks up in just 24.8 ms, whereas for PAAM at 500 ppm in over 314.0 ms.The relaxation time λ t was calculated from the filament thinning over time plotsshaded region in Fig. 2(a)-fitted with the following equation (Entov and Hinch, 1997): where R max is the maximum diameter of the filament, η p is the polymer viscosity, λ t is the longest relaxation time and τ is the surface tension.
The estimated values are listed in Table 1.The shear viscosity was measured using an ARES rheometer (ARES, TA Instruments-Waters LLC, USA) and the flow curves are shown in Fig. 2(b).The dynamic viscosity values reported in Table 1 are averages of the measured data points in Fig. 2(b).The 500 ppm PAAM solutions shows a slight shear thinning behaviour; however, the viscosity only decreases from 13.1 mPa.s to 12.0 mPa.s, when the shear rate increases from 6 s − 1 to 200 s − 1 .The maximum shear rate inside the 500 ppm Boger fluid droplets-estimated from the velocity measurements-is around 70.8 s − 1 in this study, hence the effect of shear thinning is expected to be negligible.The viscosity contrast λ between the dispersed and continuous phase ranged between

Table 1
The physical properties of the fluids used at 20 ℃.   0.1 and 0.2, i.e. it was always λ < 1.
The interfacial tension was measured by the pendant-drop method (Berry et al., 2015;Huang et al.), while the refractive index of the fluids using a refractometer (2WAJ, Shanghai Optical Instrument Co., China) as in previous studies (Li et al., 2021;Li et al., 2020).The latter were found to be equal to 1.414, 1.419, 1.472 for the 55 wt%, 66.5 wt% glycerol aqueous solutions and sunflower oil, respectively.The refractive index of PDMS is 1.412 resulting in a maximum relative refractive index difference of only 4.08%; hence the errors in reported PIV measurements due to refractive index mismatch can be neglected.
The total flow rate Q t was varied from 0 to 60 μl/min.The corresponding Reynolds number, Re = ρ c V t W/μ c based on the total flow velocity ranged from 0 to 0.12, implying that inertia forces are less important than viscous ones in this flow.ρ c and μ c are the density and shear viscosity of the continuous phase respectively, V t is the average velocity of the channel and W the characteristic length, in this case the width of the channel.The Re based on droplet velocity also varies from 0 to 0.18, confirming the absence of inertia in the flows under investigation.The Weissenberg number, which compares elastic to viscous stresses, is defined as Wi = λ t V t /W where λ t is the polymer solution relaxation time and V t the characteristic velocity.To account for the contribution of the solvent viscosity, an effective Wi is defined by multiplying by 1-β where β is the solvent viscosity ratio (Burshtein et al.,   2017); μ s /μ d , i.e.Wi eff = (1 − β)λ t V t /W ; its values ranged from 0 to 1.50.The elasticity number -indicating the relative importance of viscoelastic and inertia effects in the flow-is defined as Wi/Re (El eff = (1 − β)μ c λ t /ρ c W 2 ) and ranged from 0 to 16.14.

Microchannel geometry and experimental set up
A T-shaped microchannel (Fig. 3) was employed for the formation of microdroplets in the present experiments due to ease of fabrication and control (Thorsen et al., 2001); compared to cross-shaped microchannels, this geometry results in the formation of larger viscoelastic droplets at high Q t because of shearing stresses acting on one side of the filament.
The channel width was W = 300 μm, the continuous and dispersed phase inlets were 300 μm and 150 μm respectively, whereas the channel height H = 120 μm, as shown in Fig. 3.In this experiment, the dimensionless droplet size was l * d = l d /W = 1.5, where l d is the length of the droplets.The microchannels were fabricated in PDMS using standard soft lithography methods as described in our previous research (Li et al., 2020;Liu et al., 2017).
The experimental system used in our previous work on flow inside droplets (Ma et al., 2014), was utilized.The PDMS microchannel was placed on an inverted microscope platform (Leica DM ILM, Germany) and perfused using a syringe pump (Fusion 4000, Chemyx Inc., USA).
Glass syringes (2.5 ml, Hamilton1000, Switzerland) were used to minimize any compliance interference due to the high pressure inside the syringes.Nile red fluorescent particles (FluoSphere F-8819, Invitrogen Inc., USA) with a diameter of 1 μm were added to the dispersed phase at a concentration of 3.3% v/v.The flow was illuminated by a dual-pulse 532 nm laser Nd: YAG (New Wave, USA) and the particle fluorescence was captured by a CCD camera (C8484-05CP, Hamamatsu, Japan) in a dual-frame mode at 6 Hz through a 10 × objective lens.The time interval of the double frames was set to 0.4 ~ 4 ms according to Q t .The image resolution was 1344 × 1024 pixels 2 .LabVIEW was used for system control and image acquisition.The focal plane was set at channel half height, where fluorescent particles at the front and back ends of the droplet could be clearly displayed, as shown in Fig. 3. 800 to 1000 images were acquired, and more than 20 images containing whole droplets were selected for the calculation of average flow fields; the latter was obtained from 40 droplets whose size varied up to 2.0% for the maximum polymer concentration (material no.C) and highest Q t (45 μl/ min) employed.
Custom made Matlab routines were used to shift the droplets at different positions in the image to the same center of mass (Liu et al., 2017).Dynamic Studio (DANTEC Dynamics) was used to calculate the velocity vectors, and processing steps included masking and adaptive cross correlation with a starting interrogation window of 64 × pixels 2 , final one 32 × 32 pixels 2 and a step size of 16 × 16 pixels 2 .Vector validation was performed using an average filter and moving average with a window of 3 × 3 pixels 2 .The estimated absolute velocity fields inside the droplets were converted to the frame of reference of the droplet by subtracting the droplet velocity V d .

Effect of Wi eff on internal flow characteristics
The flow rate Q t determines the magnitude of viscous shear stresses and pressure drop in the microchannel flow and plays a key role on droplet shape and internal flow topologies (Jakiela et al., 2012;Li et al., 2021;Li et al., 2020;Kinoshita et al., 2007;Hodges et al., 2004;Liu et al., 2017;Ma et al., 2014;Kovalev et al., 2018).In viscoelastic fluids, elastic forces arising from polymer stretching, come into play and their magnitude is also affected by Q t .To study the effects of the competitive balance between elastic and viscous forces on the flow topology inside Boger fluid droplets, a 500 ppm PAAM polymer solution was used (fluid C in Table 1) and the Wi eff was controlled by adjusting Q t .
Fig. 4(a) to (c) show the streamlines inside Boger fluid droplets with increasing Wi eff .In order to identify the vortical regions within the droplet we apply the Q-criterion proposed by Hunt (Hunt et al., 1988).Q is the strain rate tensor and Ω = 1 2 [∇U − (∇U) T ] the rotation rate tensor.Positive Q regions correspond to rotation dominated regions whereas negative ones to straining regions of the flow.Fig. 4(d) to (f) show contours of the positive Q values normalized by their maximum values, i.e.Q * = Q/Q max , superimposed on the measured vector plots.
At Wi eff = 0.04, there are four eddies inside the droplets; the eddies on the upper part of the droplets rotate clockwise, whereas the lower part ones in the opposite direction.The droplet boundary velocity is consistent with the droplet movement direction, as the vectors in Fig.(d) indicate.When Wi eff increases to 0.17, the eddies at the front of the droplet expand whereas the ones at the rear shrink.There is no change in the direction of rotation, as shown in Fig. 4(e); however, the velocity direction near the droplet wall boundary reverses, leading to a new pair of eddies on either side.When Wi eff further increases to 0.33, the eddies in the front of the droplet become so strong that the back eddies disappear.The flow topology in Fig. 4(b) indicates an intermediate state as the flow transitions from a four eddy to a two-eddy structure.It is worth noting that these flow topologies have also been found in Newtonian and shear thinning droplets; the same order of flow transitions has been observed therein, especially inside Newtonian droplets, where a similar intermediate transition state has been observed (Li et al., 2020;Gu and Liow, 2011).Jakiela et al. (Jakiela et al., 2012) found that additional secondary eddies make their appearance at the caps of the droplets as Ca increases resulting in a six eddy flow structure; the differences in the number of eddies and the order of transitions can be attributed to the aspect ratio of the channel and the Ca range studied; (Li et al., 2020) the low viscosity of the continuous phase selected by Jakiela et al. (Jakiela et al., 2012) allows a lower starting range of Ca, while the higher aspect ratio compared to the present work results in a stronger wall resistance in the focus plane (Kinoshita et al., 2007).According to Mießner et al. (Mießner et al., 2020) the wall resistance determines the main vortices inside the droplets whereas the additional ones are due to the gutter flows.In order to probe further the role of elastic forces on the flow topology, higher Wi eff were investigated.Fig. 6 shows the flow transitions inside the Boger fluid droplets with Wi eff increasing from 0.67 to 1.50.At Wi eff = 0.67, the flow streamlines indicate the same topology as in Fig. 4c (i.e.Wi eff = 0.33), dominated by the appearance of two main eddies, shifted towards the front of the droplet (Fig. 6(a)).However, the small regions of positive Q * near the channel wall (indicated by the red rectangle in Fig. 4(f)) disappear in Fig. 6(d) due to a thicker oil film.When Wi eff increases to 1.00, the flow topology changes significantly; two stagnation points appear on the droplet centreline (green dots in  extremely high viscosity ratio μ c /μ d ≈ 968.7 system and was attributed to the influence of the strong shear forces exerted by the continuous phase.In the present case, the viscosity contrast is much lower and remains below 1 while Wi eff approaches 1, which implies that elastic forces become comparable to the viscous ones and should come into effect.As Wi eff increases further to 1.50, the flow at the front of the droplet reverses forming two recirculation regions on either side of the centreline, as shown in Fig. 6(c) and (f).Fig. 6(d) to (f) show that with increasing Wi eff , the regions of strong rotation in the front of the droplet shift towards the centreline and cancel out, generating initially a stagnation and subsequently a recirculating flow region.As Wi eff exceeds 1, we would expect elastic forces to play a role in the flow transitions observed in Fig. 6.
A thicker oil film implies higher boundary velocity (Li et al., 2020;Liu et al., 2017), especially at the front of droplets, in the horizontal direction.As shown in Fig. 7(a), the relative axial velocity V x near the channel wall decays inside the droplets, and becomes negative under the action of viscous dissipation, wall confinement and fluid mass continuity.A clockwise vortex appears in the upper part of the droplets as a result.As Q t increases further, the absolute axial velocity near the wall becomes so high (21.6 mm/s) that after subtracting V d , V x cannot reverse and become negative, as the green arrows at the front of the droplet indicate in Fig. 7 A flow topology similar to that in Fig. 6(c) has been found inside shear thinning droplets; the regions of strong shear and high velocity were located in the rear part of the droplets due to shear thinning caused by the continuous phase (Li et al., 2021), contrary to the Boger fluid topology shown in Fig. 7(c).The viscoelastic properties of the shear thinning solutions were not fully characterized and the contribution of the elastic forces to the observed flow topology could not be determined.In order to explore the origin of the flow transition depicted in Fig. 6(c) and the role of elasticity, the droplet flow topology in the absence of polymers (i.e.non-elastic Newtonian droplets, fluid A) was studied.droplets at Ca ranging from 2.2 × 10 -2 to 1.33 × 10 -1 .The shear viscosity and the viscosity contrast of these droplets are almost half of those for the Boger fluids.The four-eddy structure observed with Boger droplets at low Q t (Fig. 4a) cannot be seen as the Ca range shown is one order of magnitude larger here.Fig. 8(a) to (c), corresponding to Ca range from 2.2 × 10 -2 to 4.4 × 10 -2 , illustrate a transition similar to that with Boger fluid droplets, i.e. a gradual disappearance of the two main eddies and the emergence of two stagnation points and a low velocity region therein.The latter gradually expands as Ca further increases to 9.9 × 10 -2 (Fig. 8(c) to (e)) with the two stagnation points (noted with green in Fig. 6(b)) moving further apart.The flow topology of the inelastic droplets seems to stabilize at the highest Ca (Fig. 8(f)) and not transition further to the one observed for the purely elastic droplets.
The differences in the flow structure between inelastic and purely elastic droplets are also illustrated in Fig. 9(a) and (b) plotting the normalized axial velocity (V x /V d ) profiles along the droplet centreline, with the x-origin located at the rear end of the droplets.The double eddy structure is accompanied by fairly constant negative velocities along the centreline.For the inelastic fluids this occurs only for the lowest Ca shown but persists for higher flow rates for the elastic ones.The flow transition is marked by the centreline velocity increasing sharply and crossing the zero axis (twice).This happens at lower Ca for the inelastic droplets for which both the peak positive velocities V x /V d and the distance between the zero crossing points increase substantially with Ca (Fig. 9(b)).On the contrary, for the elastic fluid droplets the Ca change from 6.6 × 10 -2 to 9.9 × 10 -2 , is accompanied by a decrease in the peak V x /V d value (blue and green lines, Fig. 9(a)).The differences in the centreline profiles can thus be attributed to both wall confinement effects and fluid viscoelasticity.
Fig. 10 shows that the mean vorticity magnitude inside the droplets increases almost linearly with Q t for both elastic and inelastic droplets.However, once Ca > 9.9× 10 -2 , the vorticity inside 500 ppm Boger fluid droplets exceeds that of the inelastic droplets (see purple dotted ellipse), as a result of the observed transition.It should be noted that even if Ca is increased to 1.33 × 10 -1 , the flow topology inside the Newtonian droplets (Fig. 8(f)) does not transition to that inside the Boger droplets displayed in Fig. 6(c).Unfortunately the Ca of 1.33 × 10 -1 could not be reached with the 500 ppm PAAM solutions as Boger droplets with l * d = 1.5 could not form stably. Based on the above comparisons it can be inferred that the differences in flow topology between Figs. 6(c) and 8(e) can be attributed to the presence of the PAAM polymers.

The role of El eff on the flow topology inside the Boger droplets
In order to further probe the role of elastic forces on the measured flow topologies -hence explain the differences between Figs. 6(c) and 8 (e)-we varied El eff which is a property of the fluids only and does not depend on the shear rate.We did this by varying the concentration of PAAM in the polymer solutions; fluid B in Table 1 corresponding to a lower concentration of PAAM-and hence lower relaxation time-was selected to generate droplets at the same Q t (45 μl/min, Ca = 9.9 × 10 -2 ) as in Figs.6(c) and 8(e).This allowed to relate changes in the flow fields inside droplets to El eff .The shear viscosity of the lower PAAM concentration fluid is closer to the inelastic Newtonian one.Both Boger fluids exhibit El eff numbers higher than 1 implying that elasticity effects should dominate over the viscous ones.(e) (El eff = 0).The low velocity region at the front part of the droplet (as indicated by the yellow rectangle, Fig. 11(c) stretches along the centreline and towards the interior.With increasing El eff further it extends towards the interior and two recirculation regions appear on either side of the centreline, as was shown in Fig. 6(c).
The flow type parameter ξ, defined as (Zilz et al., 2012) |D|+|Ω| is commonly employed to visualize the regions of extension, rotation and shear in viscoelastic flows.Values of ξ = − 1, 0 and 1 imply that the flow is rotational, shear or extensional motion dominated, respectively.Fig. 12(a) to (c) show contours of ξ for increasing El eff , i.e. corresponding to the streamlines shown in Fig. 8(e), 11(a) and 6(c) respectively.The Newtonian flow exhibits a significant region of extensional flow in the middle part of the droplet surrounded by regions of rotational flow concentrated mainly at the rear side boundaries of the droplet.As El eff increases, the extensional flow at the front of the droplet is replaced by regions of shear flow as the fast moving flow around the droplet boundary meets the slow-moving flow along the centreline (black ellipse, Fig. 12(b)), which corresponds to Fig. 11(b).A further rise in El eff results in the creation of two new regions of rotational flow at the front region of the droplets corresponding to the pattern identified in Fig. 6  (c).
More quantitative comparisons of ξ can be made through selected ξ profiles plotted along the dashed line shown in Fig. 12(c).Fig. 13 illustrates the changes ξ undergoes with increasing El eff .The large region of extensional flow observed for Newtonian droplets shrinks with El eff concentrating around a narrow region around the centreline at the front of the droplet -where most of the topological changes occur-with rotational flow becoming more dominant on either side of the centreline (blue line) for El eff of 16.14, as shown in Fig. 13.

Decoupling viscosity contrast and elastic effects on the flow fields
Although elasticity is the main difference between the Boger and Newtonian fluids employed here, elasticity and viscosity vary simultaneously with increasing polymer concentration; this results in a change in the viscosity contrast between continuous and dispersed phases which is known to have an impact on the flow topology (Li et al., 2020;Ma et al., 2014).To separate the effects of elasticity and viscosity on the   1) was utilised as the dispersed phase.Thus, comparing the flow inside droplets generated using 66.5 wt % glycerol/water mixtures with the 500 ppm polymer ones (fluids C and D) the effect of elasticity can be decoupled from that of the viscosity contrast.Likewise, by comparing 55 wt% and 66.5 wt% glycerol/water droplets (fluid A with D), the effect of the viscosity contrast can be examined.

Elasticity, flow transition and instabilities
The flow transition for the high elasticity case can be attributed to the well-known coupled effects of polymer stretching along curvilinear streamlines, which result in a hoop stress drawing the fluid inwards (Pakdel and McKinley, 1996).More specifically, when a polymer molecule is stretched along curved streamlines it creates a volume force normal to the direction of the flow, towards the center of the curvature, which is termed hoop stress.Fig. 6 showed that as Wi increases the flow structure transitioned from the two-eddy structure to one characterised by a region of closed, curved streamlines.Stretching of the polymers along these generates a radial component as can be seen in Fig. 15; this might amplify the hoop stresses -pointing to the inside the curved interface as the arrows show in Fig. 15(b)-and generate an axial velocity component towards the interior.
A closer look at the measured velocity fields implies that the observed flow transition for the highly elastic droplets might be associated with the onset of an elastic instability.Figs.16 and 17 plot the instantaneous velocity distribution across the droplet at x/L = 0.77 (see purple line in Fig. 12(c)) both as profiles and spatiotemporal maps.The instantaneous velocity distributions for the two inelastic cases are fairly consistent and fluctuate little around the mean.On the contrary the velocities inside the 500 ppm Boger droplets exhibit a periodic behaviour indicating that the flow becomes unstable when Wi eff increases to 1.50.Peak velocities reach higher magnitude compared to the inelastic cases and alternate on either side of the centreline.This periodicity is not associated with a change in the shape of the droplets as the droplet interface appears stable and size varied up to 2 % in the droplets selected for flow analysis.Fully characterizing this instability is not possible due to the temporal resolution of the microPIV measurements in the present study; further work with high speed microPIV imaging is required to resolve such instabilities.
Elastic instabilities arising from elastic stresses are well known in viscoelastic flows.contours of the time averaged values of M inside 500 ppm PAAM droplets with 55% glycerol/water for Ca = 2.2 × 10 -2 , 6.6 × 10 -2 , and 9.9 × 10 -2 (the corresponding Q t values are 10, 30 and 45 μl/min respectively, and the corresponding Wi eff values are 0.33, 1.00 and 1.50) and 250 ppm droplets for Ca = 9.9 × 10 -2 (the corresponding Q t and Wi eff values are 45 μl/min and 0.40).It can be seen that the M values increase significantly, above the critical value, inside the 500 ppm PAAM droplets when Ca = 9.9 × 10 -2 , reaching a maximum value of about 6.4, especially in the front part of the droplet where the flow transition appears.The critical M value for the onset of elastic instability is estimated around M = 3 in the current work.It should be noted that exploring higher fluid elasticity regimes (hence higher M) is challenging due to instabilities in droplet formation which compromises droplet size uniformity and makes the average flow field difficult to calculate.The distributions of ξ (Fig. 12) and M (Fig. 18) imply that some analogies can   be drawn between the onset of instability in the present study and those in other extension dominated or mixed kinematic set ups, for example contraction/expansion, pore and cross-slot geometries, whereby extensional stretching of polymer chains and streamline curvature alter the flow topology/local kinematics and induce instabilities (Ekanem et al., 2020;Ekanem et al., 2022;Haward et al., 2016;Xi and Graham, 2009;Davoodi et al., 2019;Carlson et al., 2022).Interfacial tension is also expected to play a role on the onset of elastic instabilities in liquid-liquid flows, as observed by Davoodi et al. (Davoodi et al., 2021) and merits further investigation.
The elasticity driven flow transition observed herein has not been previously reported in the context of microdroplets in confined microchannels-to the best of our knowledge.However, studies with Boger droplets falling through stationary Newtonian fluids have shown that elastic forces can overcome surface tension forces, deform the droplet inducing a dimple in the rear and entrain fluid into the interior for large droplets (Sostarecz and Belmonte, 2003).This entrainment phenomena can lead to a filament which can become unstable or lead to the formation of a toroidal droplet; they have been attributed to a circular stagnation line appearing at the rear of the droplet and found dependent on the viscoelastic properties and the mobility/ size of the droplets (Emamian et al., 2022).Such phenomena can be exploited to form polymeric particles with well controlled toroidal or spiral internal structures (Sharma et al., 2012) highlighting the importance of elastic flow transitions for advanced manufacturing applications; elasticity driven transitions can either occur due to the nature of the fluids employed or engineered to enhance mixing, chemical reaction synthesis or encapsulation processes.

Conclusion
The internal flow characteristics of Boger fluid microdroplets were explored and compared with those of Newtonian (inelastic) droplets by means of micro-PIV measurements.At low Wi or Ca the flow topologies inside the Boger fluid droplets are similar to those of Newtonian droplets, transitioning from a four to a double eddy structure and subsequently to a pattern characterized by the presence of two stagnation points along the centreline and no vortical structures.As the Wi increased the elastic forces become more important inducing a transition in the flow topology inside Boger fluid droplets, leading to the formation of two recirculation regions on either side of the centreline in the front part of the droplets.This flow transition was shown to be elasticity driven; it was not observed for inelastic droplets with same viscosity contrast and Ca and occurred for elasticity numbers higher than one.It was also found to be accompanied by an elastic instability which manifests with large periodic fluctuations of the measured velocity field.
The flow changes observed in elastic (Boger) microdroplets can be attributed to the hoop stresses resulting from the coupling of elastic effects and curved streamlines.Many applications involve viscoelastic fluids, being able to tune the internal structure of the droplets through viscoelasticity providing a route to control and fine tune applications such as micromixers and biochemical reactions as well to manufacture novel formulations.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Selected images of necking behavior of PAAM solution filaments of different concentrations.The initial time is set when the filament diameter suddenly changes.The frame rate of the camera is 121 Hz.

Fig. 2 .
Fig. 2. (a) Evolution of the normalised filament diameter.The inset figure is a magnified view of the filament necking process.(b) Steady shear viscosity of the dispersed phase fluids.

Fig. 3 .
Fig. 3. Schematic of the experimental setup showing details of the microchannels and the imaging system.

Fig. 5
illustrates schematically the vortex formation mechanism in the above flow topologies.According to the direction of the vectors inside the droplets, shown in Fig. 4(d) to (f), eddies O 1 and O 3 (similarly eddies O 2 and O 4 ) are caused by the action of the continuous phase shear stresses F 1 and F 2 , as shown in Fig. 5(a).As Wi eff or Ca increases, the eddies O 3 and O 4 become stronger, as shown in Fig. 5(b), implying that F 2 become larger and the region of strong rotation inside the droplets moves forward.Since the oil film is thinner towards the back of droplets compared to that in the front part as shown in Fig. 5(b), the shear stress F 3 caused by the channel wall leads to the formation of eddy O 5 and the transitional flow topology observed at Wi eff = 0.17.The rotation directions of O 1 and O 3 are the same, therefore, when F 2 (Wi or Ca) becomes large enough, they merge into one, as shown in Fig. 5(c).In addition, due to the thicker oil film at high Q t , the region of positive Q * near the upper and lower channel walls shown in Fig. 4(e) becomes smaller, as the red rectangle indicates in Fig. 4(f).Since Wi eff is always less than 1 in the range investigated (Wi eff = 0.04 0.33), the viscous forces are expected to dominate the flow topology changes and the formation mechanisms appear similar to Newtonian droplets.

Fig. 4 .
Fig. 4. (a) to (c) show the streamlines inside 500 ppm PAAM Boger fluid droplets with 55 wt% glycerol/water added atWi eff = 0.04, 0.17, and 0.33 respectively, the corresponding Ca is 2.7 × 10 -3 , 1.1 × 10 -2 and 2.2 × 10 -2 , (d) to (f) show relative velocity vectors (based on Q * contours, positive Q * regions correspond to rotation dominated regions whereas negative ones to straining regions of the flow).V d in Figs.(d) to (f) are 0.58, 2.8 and 5.6 mm/s respectively, and used as reference vectors to indicate velocity magnitude.The corresponding Q t values in Fig. 4(a) to 4(c) are 1.2, 5, and 10 μl/min respectively.The regions of positive Q * near the channel wall indicated by the red rectangle in Fig. 4(f) shrink preliminary due to a thicker oil film.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6
Fig. 6(b)).The two main vortices are no longer visible and the flow pattern is dominated by a region of low velocities and closed streamlines-directly connecting the middle and front parts of the droplets (red rectangle in Fig. 6(b)-and a jet like flow emanating from it.A similar flow topology was observed by Kovalev (Li et al., 2021) et al for an (b).Because the rotational directions of O 3 and O 4 are opposite as shown in Fig.5(c), they come close due to higher F 2, and then cancel.As a result, a region of relatively low, almost stagnated velocities around the centreline emerges, dividing the droplet into two regions of strong negative velocities at the front and rear of the droplets.Due to fluid mass continuity and wall confinement, especially for the channel wall parallel to the focus plane(Kinoshita et al., 2007), strong negative V x values in the rear part of the droplet exist resulting in the flow topology of Fig.6(b).As Wi eff increases, the high velocity flow near the wall boundaries moves towards the front of the droplet compressing the negative velocity region and forcing the flow to recirculate (Fig.6(c)).
Fig. 6.(a) to (c) show the streamlines inside 500 ppm PAAM droplets for Wi eff = 0.67, 1.00, and 1.50 respectively, the corresponding Ca is 4.4 × 10 -2 , 6.6 × 10 -2 and 9.9 × 10 -2 .(d) to (f) are relative velocity vectors (with superimposed Q * contours).V d values in Fig. 6(d) to (f) are 11.6, 18.4 and 28.2 mm/s respectively.The corresponding Q t in Fig. 6(a) to 6(c) are 20, 30, and 45 μl/min.The red rectangle in Fig. 6(b) indicates the low velocity region and stagnation points indicated by green dots.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.Mean vorticity magnitude as a function of Ca for inelastic and elastic (Boger) fluid droplets.

Fig. 11 .
Fig. 11.(a), (b) and (c) show the streamlines, relative velocity vectors (superimposed on Q * contours) and V x contours inside 250 ppm Boger fluid droplets, respectively; El eff = 4.28 andQ t = 45 μl/min.(d) V x contour inside the inelastic (55 wt% glycerol/water) droplet (El eff = 0) forQ t = 45 μl/min.V x is normalized by V d .V d = 29.9mm/s in Fig. (a), (b) and (c).V d = 30.0mm/s in Fig. (d).The yellow rectangle in (c) indicates the extension of the low (negative) velocity region along the droplet centreline, shown by purple dashed line.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 14(a)  and (b) show the streamlines and relative velocity vectors inside 66.5 wt% glycerol aqueous droplets at Q t = 45 μl/min respectively; the corresponding capillary number is Ca = 1.1 × 10 − 1 due to slightly reduced interfacial tension.The flow patterns and ξ contour appear similar to those found in 55 wt% glycerol aqueous droplets shown in Figs.8(e) and 12(a).In addition, the normalized axial centreline velocities inside the 55 wt% and 66.5 wt% glycerol/water droplets are similar.This indicates that the higher droplet viscosity may increase the flow resistance or pressure drop in the channel, but it does not affect the flow topology and the velocity distribution.Thus, it can be inferred that the viscosity contrast has no or only a weak effect on the flow fields inside droplets in the present study.Ma et al. (Ma et al., 2014) et al. reported that when the droplet viscosity contrast increases, the flow pattern becomes more uniform.Although the viscosity contrast increases by a factor of about 2 in fluid D it remains low and well below 1. Comparing Fig.14(a) with Fig. 6(c) (fluids D with C) shows differences in topology which can be attributed only to fluid elasticity.Based on the above comparisons, we can verify that the flow transitions shown in Fig. 6(b) and (c) are driven mainly by elastic rather viscous forces.
They have been shown to create secondary flows and lead to elastic turbulence in a number of flows including microscale ones (Poole, 2019).The Pakdel-McKinley criterion introduced by McKinley et al. (McKinley et al., 1996) combining the effects of streamline curvature and tensile stress in one dimensionless parameter has been found to predict the onset of purely elastic instabilities well in a number of flows, including two phase flows as was recently demonstrated by Davoodi et al. (Davoodi et al., 2021) The Pakdel-McKinley criterion is expressed as: M = ̅̅̅̅̅̅̅̅̅̅̅̅̅ τ 11 λ t v ′ d Rμ d γ √ where τ 11 is the elastic normal stress in the streamwise direction, expressed as τ 11 = 2(μ d − μ s )λ t γ2 (Casanellas et al., 2016); v ′ d is the local axial velocity, R is the local curvature of the streamlines and γ2 is the local shear rate estimate from the measured velocity field.Fig. 18 shows

Fig. 12 .Fig. 13 .
Fig. 12. (a), (b) and (c) show the ξ contours inside droplets forEl eff = 0, 4.28, and 16.14 respectively.Q t = 45 μl/min in all cases.The black ellipse in (b) marks the region where the kinematics shifts from extensional to shear, and the purple dotted line in (c) is where profiles of the flow parameter are extracted for more detailed comparison in Fig. 13.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 16 .Fig. 17 .
Fig. 16.The instantaneous velocity distribution along the green line in Fig. 12(c) for 40 droplets at 45 μl/min.(a) 55 wt% glycerol/water droplet.(b) 66.5 wt% glycerol/water droplet.(c) 500 ppm PAAM droplet with 55% glycerol added.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)