1 Introduction

Microscopic two-phase flows find application in thermal (Magnini et al. 2013), chemical (Tanimu et al. 2017), biological (Mazutis et al. 2013) and medical (Kang et al. 2014) processes. An overview of the application scope of Taylor flows is given by Chou et al. (2015). Taylor flows are often realized in rectangular flow structures, which are integrated in microfluidic chips, that feature pressure-driven flows.

Downsizing multiphase flows increase the specific surface area of the flow, promote heat and mass transfer (Bandara et al. 2015), enable precise handling of small sample volumes (Garstecki et al. 2006; Whitesides 2006) and allow high-speed processing without cross-contamination (Seemann et al. 2011). Cross-contamination is avoided by a film of continuous phase that lubricates the wall and separates the disperse phase from the wall. In rectangular channels, the thin lubricating wall films are accompanied by the so-called gutters (van Steijn et al. 2009), which enable a bypass flow through the continuous phase-filled corners (Kreutzer et al. 2005).

A good comprehension of the three-dimensional (3D) velocity distribution and its derived quantities, e.g., the shear rate magnitude and distribution, is important for optimal process design. The interaction of advection and diffusion influences the efficiency of multiphase processes, while the resilience of educt or product against shear forces limits the applicable velocity gradients (Liu et al. 2017).

For low capillary numbers Ca\(_{\mathrm {c}} = \tfrac{u_{\mathrm {0,tot}} \eta _{\mathrm {c}}}{ \sigma }\), the flow inside and around Taylor droplets is dominated by surface tension forces. Herein, \(\sigma\), \(\eta _{\mathrm {c}}\) and \(u_{\mathrm {0,tot}}\) describe the interfacial tension, the dynamic viscosity of the continuous phase and the superficial velocity, respectively. The superficial velocity \(u_{\mathrm {0,tot}}=\tfrac{Q_{\mathrm {c}}+Q_{\mathrm {d}}}{A_{\mathrm {ch}}}\) is derived from the total volume flow through the area of the channel cross section \(A_{\mathrm {ch}}\), where \(Q_{\mathrm {c}}\) and \(Q_{\mathrm {d}}\) represent the volume flow of the continuous and the disperse phase, respectively. Avoiding high shear forces necessitates reduced velocity gradients in small channel dimensions. Thus, the viscous forces dominate the inertia forces Re\(_{\mathrm {c}} = \frac{\rho _{c} u_{\mathrm {0,tot}} H_{\mathrm {ch}}}{\eta _{\mathrm {c}}} \le 1\). Here, \(\rho _{\mathrm {c}}\) denotes the density of the continuous phase and \(H_{\mathrm {ch}}\) is the channel height. Buoyancy forces are negligible, since the Bond number Bo = \(\tfrac{\varDelta \rho g H^2_{\mathrm {ch}}}{\sigma } \ll\) 1. Herein, \(\varDelta \rho\) denotes the density difference between the two phases and g is the gravitational acceleration.

A number of simulations have been published to clarify the influence of the 3D velocity distribution of pressure-driven Taylor flows in rectangular microchannels on the involved processes. Hazel and Heil (2002) investigated the propagation of a semi-infinite bubble at \(0.001 \le \mathrm{Ca}_{\mathrm {c}} \le 10\) and Re\(_{\mathrm {c}} \ll 1\). Taha and Cui (2006) performed CFD modeling of the rising slug flow inside square capillaries (\(0.001 \le \mathrm{Ca}_{\mathrm {c}} \le 3, \mathrm{Re}_{\mathrm {c}} \ll 1\)). Raimondi et al. (2008) simulated the mass transfer in square microchannels for liquid–liquid slug flow (\(0.001 \le \mathrm{Ca}_{\mathrm {c}} \le 10\) and Re\(_{\mathrm {c}} \ll 1\)), and Guido and Preziosi (2010) studied the droplet deformation of Taylor flows (\(0.001 \le \mathrm{Ca}_{\mathrm {c}} \le 1, \mathrm{Re}_{\mathrm {c}} \le 50\)). Falconi et al. (2016) presented an experimental study for vertical pressure-driven Taylor flow in rectangular minichannels (Ca\(_{\mathrm {c}} = 0.1, \mathrm{Re}_{\mathrm {c}} = 7\)). Bordbar et al. (2018) reviewed the slug flow in microchannels with respect to numerical simulation and applications. Kumari et al. (2019) published a study on the flow field and heat transfer in slug flow in microchannels with respect to the influence of the bubble volume. Recently, an analytical solution of the Taylor flow in an elongated ellipsoid droplet in the frame of 2D Stokes equations was published by Makeev et al. (2019). Their solution delivers a benchmark to evaluate numerical algorithms that aim to calculate the flow in and around Taylor droplets. The work clearly links the driving velocity distribution on the droplet interface to the circulation of the flow inside: A single vortex develops in the ‘bypass’ mode, and a set of three vortices result from the ‘slug flow’ mode.

On the other hand, detailed measurements of the velocity field in and around microscopic Taylor droplets in rectangular microchannels are reported by Kinoshita et al. (2007). They used a high-speed confocal scanning procedure that enables multiple measurement planes with 2D velocity information. An integration of the continuity equation results in the out-of-plane component. The scanning process of the high-speed confocal scanner with Nipkow disk allows measurements at the following approximated dimensionless numbers: Ca\(_{\mathrm {c}}\) = 0.0015, Re\(_{\mathrm {c}}\) = 0.00025, \(A_{\mathrm {r}}\) = 1.72, MD\(^*\) = 0.032. Herein, \(A_{\mathrm {r}}\) denotes the aspect ratio relating the channel width to its height, while MD\(^*=\tfrac{\text {MD}}{H_{\mathrm {ch}}}\) refers to the measurement depth MD of a single optical slice of the experimental setup divided by the channel height \(H_{\mathrm {ch}}\). The same research group extended their technique and reported a simultaneous measurement of a Taylor droplet including its surroundings (Oishi et al. 2011). An approach to distinguish between tracer particles of different fluorescent color enabled them to discriminate between the phases. The measured Taylor flow features dimensionless numbers, which are comparable to those of our study (Ca\(_{\mathrm {c}}\) = 0.003, Re\(_{\mathrm {c}}\) = 0.001, \(A_{\mathrm {r}}\) = 2.44, MD\(^*\) = 0.084). Mohammadi and Sharp (2013) summarized the experimental methods to access the velocity distribution of Taylor flows. Meyer et al. (2014) investigated vertically rising Taylor bubbles in a rectangular mini channel. They measured the velocity field of the continuous phase in the slugs and in the film to support a numerical simulation that provides the entire flow. (Ca\(_{\mathrm {c}}\) = 0.1, Re\(_{\mathrm {c}}\) = 8, \(A_{\mathrm {r}}\) = 1, MD\(^*\) = 0.015). Khodaparast et al. (2013) developed a measurement system (microparticle shadow velocimetry—µPSV) to investigate Taylor flow in circular microchannels in detail (Khodaparast et al. 2014). The dimensionless parameters are: Ca\(_{\mathrm {c}}\) = 0.0003–0.09, Re\(_{\mathrm {c}}\) = 0.1–10, MD\(^*\) = 0.071. The measurements of Ma et al. (2014) revealed a viscosity ratio-dependent change of flow topology inside their droplets that move in rectangular microchannels. They base their findings on the flow field measurements of two focal plane positions (Ca\(_{\mathrm {c}}\) = 0.001–0.1, Re\(_{\mathrm {c}}\) = 0.1–10, \(A_{\mathrm {r}}\) = 1.33, MD\(^*\) = 0.099). Liu et al. (2017) provide a µPIV investigation of the internal flow transitions inside droplets traveling in a rectangular microchannel. Their measurements sample the 3D velocity field solely in the symmetry plane (Bo = 0.0043, Ca\(_{\mathrm {c}}\) = 0.005–0.26, Re\(_{\mathrm {c}}\) = 0.1–10, \(A_{\mathrm {r}}\) = 1.5, MD\(^*\) = 0.090).

In the present study, we experimentally capture the 3D flow field of a Taylor flow at low Ca\(_{\mathrm {c}}\) and Re\(_{\mathrm {c}}\). The 3D velocity data may serve as benchmark for the development of new theoretical and numerical models for Taylor flow and for parametric CFD studies. The hydrodynamic features of Taylor droplets are investigated in a horizontal square microchannel (Ca\(_{\mathrm {c}}\) = 0.005, Re\(_{\mathrm {c}}\) = 0.051, \(A_{\mathrm {r}} \approx 1\), MD\(^*\) = 0.056). Based on a multi-plane µPIV measurement, we simultaneously measure the 3D2C velocity field inside and outside of Taylor droplets. Subsequently, we reconstruct the full 3D3C velocity field using the continuity equation. We show the stagnation regions of the flow relative to the droplet motion and deduce the normalized shear rate distribution. This delivers information about  the location of  maximum strain. The position is of interest for handling delicate products or sensitive educts. Our novel approach yields an experimental data set that allows the visualization of the streamlines of the velocity distribution on the interface of a Taylor droplet and to directly relate it to the main and secondary vortices of the velocity field. We provide the measurement data as a digital resource for the validation of numerical simulation or further assessment.

2 Materials and methods

At first, we introduce the materials and methods used to measure the 3D2C velocity distribution of a Taylor flow by means of µPIV. Secondly, we reconstruct the third velocity component and verify the 3D flow distribution. Subsequently, we present the methods for the shear rate quantification and the vortex visualization.

2.1 Experimental setup

Syringe pumps provide stable and pulsation-free volumetric flow rates of continuous phase \(Q_{\mathrm {c}}\) (octanol) and disperse phase \(Q_{\mathrm {d}}\) (water–glycerol) to a flow focusing (FF) junction molded in polydimethylsiloxane (PDMS) (Mießner et al. 2008). The flow rate ratio is \(\alpha = Q_{\mathrm {d}}/Q_{\mathrm {c}} = 1\) and delivers a disperse phase fraction of \(\varepsilon _{\mathrm {d}} = \tfrac{Q_{\mathrm {d}}}{Q_{\mathrm {d}}+Q_{\mathrm {c}}} = \tfrac{\alpha }{\alpha +1} = 0.5\). A regular Taylor droplet train forms behind a pinhole. The droplets are always separated from the wall by a thin octanol film of the thickness \(\delta \lesssim\)\(1\, \upmu \hbox {m}\).

The cross section of the microchannel has a slight trapezoidal deviation from the square shape (\(W_{\mathrm {ch,1}} = {104}\, \upmu \hbox {m}, W_{\mathrm {ch,2}} = {96}\, \upmu \hbox {m}, H_{\mathrm {ch}} = {96}\, \upmu \hbox {m}\)). The cross section forms a convex isosceles trapezoid with channel sidewalls enclosing an angle of \({4.7}^{\circ }\) in comparison with an angle of \({0}^{\circ }\) between the parallel sidewalls of a square microchannel. The trapezoidal deviation from an ideal square aspect ratio \(A_{\mathrm {r}} = 1\) is a result of the fabrication process. We consider this shape deviation to be small and treat the cross section as rectangular.

SEM imaging and visual inspection suggest that the surface roughness (\(\mathrm {R}_{\mathrm {A}} \lessapprox {0.2}\, \upmu \hbox {m}\)) is small in relation to the channel height. The disperse phase consists of an aqueous mixture of water and glycerol. Changing the mass fraction to about 18% water and 82% glycerol allows the adjustment of the refractive index (RI) of the disperse phase to that of the continuous phase (octanol, RI=1.428). A more detailed description of the refractive index matching is given by Mießner et al. (2019). The flow in and around the Taylor droplets is optically measured 5 mm downstream of the FF junction.

The continuous phase as well as the droplets are seeded with fluorescent tracer particles. The continuous-phase seeding consists of polystyrene (PS) (microparticles GmbH) with a nominal particle diameter of \(d_{\mathrm {p}} ={0.9}\, \upmu \hbox {m}\), while the aqueous disperse phase is seeded with polyethylenglycol (PEG)-coated PS particles (\(d_{\mathrm {p}} = {1.3}\, \upmu \hbox {m}\)). The seeding density of the continuous phase is set higher than the seeding density of the droplet to enable discrimination between both phases. All particles are coated with Rhodamine B, which emits a fluorescent light (nominal wavelength at maximum intensity: \(\lambda _{\mathrm {peak}} = {600}\, \hbox {nm}\)) upon excitation with a Nd:YAG laser pulse (\(\lambda _{\mathrm {ex}} = {532}\, \hbox {nm}\)). The double-pulsed laser light is coupled into the light path of the inverted microscope (Zeiss Axiovert 100) and a dichroic mirror that passes only the fluorescent signal to the double-frame CCD camera (PCO, Sensicam QE). Camera and laser are synchronized with the programmable timing unit of the PIV system (LaVision). A piezo-stepper accurately controls the z-position of the objective (Zeiss LD LCI Plan-Apochromat 25x/0.8) and its focal plane. We use this setup to scan through the flow volume in the z-direction with a measurement depth (MD) of about \({5.4}\, \upmu \hbox {m}\). A more detailed description of the experimental setup, the two-phase system and the refractive index matching is given by Mießner et al. (2019). A summary of the flow conditions is given in Table 1.

Table 1 Overview of parameters that influence the flow conditions of the Taylor flow
Full size table

2.2 Image acquisition and image processing

The symmetry of the Taylor flow in straight rectangular microchannels allows to conduct measurements in one half of the channel to reduce the amount of processed data. After initiating the flow, it takes about 2 minutes for droplet size fluctuation to reduce to a reasonable level. We decided to apply a tenfold time to make sure the initial fluctuation is gone. After this transient of 20 min, we consider the flow to be quasi-stationary: For the moving droplets, the flow inside and around each droplet is similar for all droplets. The µPIV setup is not synchronized to the flow and freely records a number of 600 double frames per measurement plane. The measured Taylor droplet train exhibits a droplet size distribution that varies 10% in length. The recording frequency of PIV double frame images \(f_{\mathrm {rec}}\) is fixed at 12 Hz, while the droplet frequency oscillates around \(f_{\mathrm {d}} \approx\) 14 Hz. The temporal distance between the double frames is set to 1.75 ms. The depth of correlation determines the measurement depth of the optical setup to \({5.4}\, \upmu \hbox {m}\). We divide half of the channel height into 11 equidistant (\(\varDelta H = {5.4}\, \upmu \hbox {m}\)) measurement planes to avoid recording redundant flow information. The chosen z-positions of the focal planes include the center symmetry plane and a measurement plane directly at the bottom cover of the microchannel. Beginning with the center plane, the µPIV system is programmed such that after recoding 600 images, the piezo-stepper moves the objective and thus its focus to the next plane. The procedure is repeated until the measurement is completed at the bottom lid position of the microchannel. Figure 1 shows a sketch of the measurement plane position with respect to the Taylor droplet interface.

Fig. 1
figure 1

Sketch of the measurement planes on which ensemble correlation µPIV is performed. For clarity, only 6 out of 11 measurement planes are depicted. The approximate Taylor droplet model (Mießner et al. 2019) indicates the position of the interface (light blue surface). The Taylor droplets translate in the positive x-direction, and the droplet’s front tip is placed at the origin of the Cartesian coordinate system (yellow dot)

Full size image

Due to the asynchronous recording and the strong length variability of the observed flow, preprocessing of the recorded data is applied to receive a quasi-stationary representation of the velocity distribution inside and around an average Taylor droplet. The preprocessing involves image selection (e.g., a statistical approach to reduce the droplet length variation), processing (in-frame image shifting in the x–y-direction) and vertical alignment (in the z-direction). The method is explained in detail in Mießner et al. (2019). The preprocessing reduces the droplet length deviation from 10% to about 3%. This increased precision enables the determination of the droplet velocity at the stagnation points and allows to accurately change to a moving frame of reference with respect to the droplet motion.

For each measurement plane, the identified double images of suitable droplets are selected such that at least 60 image pairs are further processed. To increase the data density, we make use of the flow symmetry along the center axis of the measured Taylor droplets. The selected image set is extended by additionally processing the mirrored data set in the µPIV analysis. This results in at least 120 post-processed image pairs per measurement plane, which are subsequently evaluated with an ensemble correlation PIV algorithm.

2.3 Ensemble correlation µPIV

On each measurement plane, the selected and processed image pairs are evaluated using ensemble correlation (Delnoij et al. 1999). The method relies on averaging coinciding correlation planes from a sequence of image pairs of a stationary flow (Raffel et al. 2018). After passing through the post-processing, the images fulfill this requirement. Stochastic influences such as Brownian motion are suppressed since randomly deviated correlation peaks are averaged out in favor of the recurring correlation peaks of the stationary flow. Additionally, this technique increases the information density of the sparsely seeded flows (Adrian and Westerweel 2011). Depending on the particle density, Vennemann et al. (2006) estimated 50 image pairs to be sufficient to obtain close to 100% reliable vectors (for interrogation areas of \(64 \times 64\) px).

A cycle of decreasing square interrogation window sizes [256 128 96 64 32 16 16] px was used with an overlap of 50% applying ensemble correlation for each iteration. The resulting vector field of each cycle is used as input information for the next iteration step. The last evaluation step is performed twice to further increase the accuracy of the resulting ensemble correlated vector field.

The depth of correlation (DOC) without shear is \(z_{\mathrm {corr,0}} = {5.4}\, \upmu \hbox {m}\) (Olsen and Adrian 2000; Lindken et al. 2009; Rossi et al. 2010). However, the flow in our setup is subject to shear: Olsen (2008, 2010) introduced an approximation to describe the shear dependency of the DOC. In-plane shear increases the DOC, while to a lesser extent out-of-plane shear decreases the DOC. Based on the statistics of the in-plane shear (Fig. 2), we found that the DOC increase in our setup ranges between 1.21\(\times\) and 2.35\(\times\) at most, depending on whether 68% or 95% of the shear rate distribution data are included in the estimation.

Fig. 2
figure 2

Deviation of the DOC due to shear. Based on the probability distribution of the in-plane shear rates (red) on the left y-axis, the cumulated probability on the right y-axis allows to estimate the amount of included values of the shear rates distribution. The bins of the shear rate distribution (bottom x-axis) are converted into the shear related depth of correlation (DOC) variation (top x-axis) (Olsen 2008), to simplify the reading of the graph

Full size image

An analysis of the correlation between the in-plane shear distribution of the downstream velocities and the reconstruction error distribution reveals that the error increases with rising shear rates. In fact, the increase could be explained with the shear increase in the DOC. In any case, given the relative low mean absolute reconstruction error presented in Sect. 2.5, we consider the change of the DOC due to shear to be acceptable and the chosen measurement plane distance to be valid.

Evaluating the particle image data sets of all measurement planes results in a 3D2C velocity distribution with an in-plane spatial resolution \(\varDelta x = \varDelta y = {1.8}\, \upmu \hbox {m}\) and an out-of-plane distance of \(\varDelta z = {5.4}\, \upmu \hbox {m}\). The number of measurement nodes of 171 \(\times\) 71 \(\times\) 11 amounts to 133551 regularly spaced vectors in one of the symmetric halves of the quasi-stationary Taylor flow.

2.4 Reconstruction of the velocity z-component w

The measurement delivers a 3D2C velocity field, which lacks the third velocity component w in the z-direction. A reconstruction retrieves the missing component based on the available measured information. We apply the continuity equation for incompressible flows to the 2D in-plane velocity components u and v. A second-order central difference scheme delivers the derivatives of the 2D velocity field (Adrian and Westerweel 2011).

Beginning with the supposedly divergence-free symmetry plane (i.e., with \(\hbox {d}w/\hbox {d}z = 0\) and \(w = 0\)) of Taylor droplets and proceeding toward the bottom cover of the microchannel, the integration of the divergence in the z-direction delivers an estimate for the out-of-plane velocity component w. A more detailed description of the method is given by Robinson and Rockwell (1993) and Brücker (1995, 1997). We do not impose a boundary condition for the velocity in the z-direction at the cover, but use the deviation from \(w = 0\) to assess the accuracy of the reconstruction.

Prior to reconstruction, an approach based on the normalized median test (Westerweel 1994; Westerweel and Scarano 2005) is applied to detect, remove and replace outliers of the measured velocities u and v and their derivatives. Additionally, a 2D median filter with a \(3 \times 3\) neighborhood for velocities and a \(5 \times 5\) neighborhood for the derivatives is applied to suppress random noise (Westerweel 1993; Adrian and Westerweel 2011).

2.5 Validation of the reconstructed 3D3C flow field

The transverse velocities v (measured) and w (reconstructed) are considered to be symmetric for a perfect aspect ratio of \(A_{\mathrm {r}} = 1\). Therefore, a comparison of the probability density distribution (PDF) of the transverse velocities (Fig. 3) offers a visual qualitative impression of the reconstruction accuracy. The shape of the probability density function representing the reconstructed velocities (red solid line) is similar to that of the measured component (black solid line). The deviations between the peaks and the flanks are attributed to the measurement uncertainties and the reconstruction.

The PDF of the normalized stream-wise relative velocity (black dotted line) represents measurement data only. The two peaks at a normalized velocity amplitude of − 1 and 0 are related to the velocities at the wall and the stagnation regions of the flow, respectively. The broad positive and negative flanks of the graph result from the bidirectional stream-wise flow of the main wall-driven vortices inside and outside the Taylor droplet.

Fig. 3
figure 3

Normalized velocities u (stream-wise direction, measurement), v (transverse direction, measurement), w (transverse direction, reconstruction). Qualitatively, the PDF of the reconstructed transverse velocity is similar to the measured transverse component PDF. The deviations at the peak and the flanks are at first attributed to the measurement uncertainties, the median filter application and the calculation of the derivative used for the reconstruction

Full size image

Secondly, the symmetry of the flow allows the comparison of the measured y-components in the first quadrant of the cross section (\(+\) y, \(+\) z) with the reconstructed z-component of the second quadrant (− y, \(+\) z): Rotating the measured 2D velocity field by \(+90^{\circ }\) and performing a linear interpolation of the rotated measured transverse velocity component on the grid of the reconstructed velocities gives \(v^{\perp }\). Subtracting the measured values from the reconstructed velocity components in quadrant two results in a 3D residue distribution, which serves as an estimate of the reconstruction accuracy:

$$\begin{aligned} {\mathcal {R}} = \frac{\left( w-v^\perp \right) }{u_{\mathrm {d}}} \end{aligned}$$

Figure 4 shows the mean residue \(\bar{{\mathcal {R}}}_z\) and its standard deviation in dependence of the channel height \(\frac{z}{H_{\mathrm {ch}}}\). The residue grows with distance to the xy-symmetry plane. The maximum deviation can be found at a channel height of \(\frac{z}{H_{\mathrm {ch}}} \approx 0.30\), i.e., the location of maximum out-of plane velocities in the z-direction. The overall reconstruction error (Fig. 4a) is estimated to be \(\bar{{\mathcal {R}}}_z \pm 1\sigma ({\mathcal {R}}_z) = 0,6 \pm 5.3\%\).

Fig. 4
figure 4

Estimation of the mean residue \(\bar{{\mathcal {R}}}_z\) distribution and its standard deviation in dependence of the channel height \(\frac{z}{H_{\mathrm {ch}}}\). Inset a) PDF of the total residue distribution in the domain \(\bar{{\mathcal {R}}}_z \pm 2\sigma ({\mathcal {R}}_z) = 0,6 \pm 10.6\%\)

Full size image

Finally, an analysis of the stream-wise and transverse volume flow inside and outside the Taylor droplets in Fig. 5 confirms the coherence of the measured flow data. The mean volume flux distribution \(Q_x, Q_y\) and \(Q_z\) along the center axis of the channel is calculated with the reconstructed 3D3C absolute velocity field and the corresponding reference areas d\(A_x = \varDelta y \varDelta z\), d\(A_y = \varDelta x \varDelta z\) and d\(A_z = \varDelta x \varDelta y\). The mean volume fluxes in the x-direction (Fig. 5a), the y-direction (Fig. 5b) and the z-direction (Fig. 5c) are non-dimensionalized with the superimposed total volume flow \((Q_{\mathrm {tot}} = Q_{\mathrm {d}} + Q_{\mathrm {c}})\). The total volume flow of a direction is marked blue, while the continuous phase and the disperse phase are shown as black and red lines, respectively. The volume flux \(Q_{i,\mathrm {c}}\) and \(Q_{i,\mathrm {d}}\) of each of the two phases is extracted by applying the approximate interface model as a logical discriminator (Mießner et al. 2019). The Taylor flow is oriented in the positive x-direction, and the Taylor droplet’s front cap tip is located at the origin of the coordinate system. The tip of the back cap is located at \(\frac{x}{H_{\mathrm {ch}}} = -1.58\), which corresponds to the dimensionless droplet length \(\frac{L_{\mathrm {d}}}{H_{\mathrm {ch}}}\).

Fig. 5
figure 5

Analysis of the dimensionless volume flow in stream-wise \(Q_{x,i}/Q_{\mathrm {tot}}\) and in the transverse direction inside and outside the Taylor droplets (\(Q_{y,i}/Q_{\mathrm {tot}}\) and \(Q_{z,i}/Q_{\mathrm {tot}}\)) to demonstrate the coherence of the flow data. Depicted are volume flows occurring through cross sections perpendicular to the main flow direction (along the x-axis). The total volume flow \((Q_{\mathrm {tot}} = Q_{\mathrm {d}} + Q_{\mathrm {c}})\) is expected to match the measurement volume flow measurement in the x-direction (\(Q_x/Q_{\mathrm {tot}} = 1\)), while the volume flows in the transverse direction are expected to be zero due to conservation of mass (\(Q_y/Q_{\mathrm {tot}} = Q_z/Q_{\mathrm {tot}} = 0\)). The indices \(_{\mathrm {d}}\) and \(_{\mathrm {c}}\) refer to the disperse and continuous phase, respectively. The stream-wise volume flux a and the transverse volume flux in the y-direction b represent the results of the 3D2C µPIV measurement, while c depicts the reconstructed transverse out-of-plane volume flux in the z-direction. The stream-wise flux deviation is attributed to the channel deformation \((94.44\% \pm 1.2\%)\,Q_{\mathrm {tot}}\), while the deviation in the y-direction highlights the accuracy of the µPIV measurement. The reconstruction of the third velocity component carries a standard deviation of ± 2.7% of \(Q_{\mathrm {tot}}\) along the longitudinal axis of the droplet

Full size image

The measured mean stream-wise volume flux distribution along the x-axis is shown in Fig. 5a. The blue graph corresponds to the superficial volume flux in the x-direction, while the black and red graphs depict the mean continuous and mean disperse phase flux per channel cross section, respectively. The mean superficial volume flow deviates with \(-5.56\%\) at most from the nominal flow rate, which is set by the syringe pumps (\(Q_{\mathrm {tot}} = {3.0}\, \upmu \hbox {l} \hbox {min}^{-1}\)). Ruling out defective equipment or leakage flow, we expect the deformation of the channel material (polymethyldisiloxane - PDMS) to be the reason for the deviation. Taylor droplets are found to exert pressure peaks to the wall while passing (Abate et al. 2012), and the ductility of the channel material allows the wall to subside, while being pressurized (Gervais et al. 2006). Thus, the pressure drop-related channel expansion causes a decrease in the mean velocity, and as a result the microchannel dilates. When recalculating the volume flux \(Q_i\) from the flow field by applying the reference areas d\(A_i\), the volume flow per cross section appears to be too small. Assuming a correct measurement procedure, this relates to an estimated cross section expansion rate of 2.9% in channel width (\({3.3}\, \upmu \hbox {m}\)) and height (\({2.8}\, \upmu \hbox {m}\)).

Proceeding along the x-axis toward the droplet’s back cap (Fig. 5a), the volume flux of the continuous phase (black graph) decreases as the droplet phase (red graph) occupies more space. The volume flux in the gutter and in the wall film is present as well: The dimensionless volume flow of continuous phase never drops below 2.1% at \(\frac{x}{H_{\mathrm {ch}}} \approx -1\). Toward the droplet tip, the volume flux development reverses until only continuous phase is present from \(\frac{x}{H_{\mathrm {ch}}} > 0\) on.

Figure 5b shows the mean transverse volume flux in the y-direction along the Taylor droplet main axis, that is expected to be zero. All fluctuations in the y-direction remain below 0.1%. This consideration reflects the measurement error introduced by the applied µPIV measurement technique.

The mean volume flow based on the reconstructed (dashed lines) transverse velocity component in the z-direction is shown in Fig. 5c. Due to the symmetry of Taylor flows, the mean flux per cross section is expected to cancel out. However, the mean residual flux amounts to \(-\,0.001\% \pm 2.7\%\) of \(\,Q_{\mathrm {tot}}\). The maximum deviation amounts to − 6% of \(Q_{\mathrm {tot}}\). The result reflects the accuracy of the applied reconstruction method, since we carefully positioned the measurement planes with respect to the symmetry planes of the microchannel. The subdivision of the channel volume into disperse and continuous phase reveals that the deviation is located in both phases.

The error sources are identified to be the data filter (low-pass filter for data noise) and the determination of the derivatives (which acts as a high-pass filter for data noise) in addition to the dilation of the PDMS microchannel and the shear influence on the DOC. From this analysis, the standard deviation of the reconstruction error is \(\pm 2.7\%\) of the total volume flux \(Q_{\mathrm {tot}}\), which is half as big as the reconstruction error estimation in Fig. 4. We expect the difference to result from the integration step taken to derive the absolute volume flow.

2.6 Shear rate distribution

Based on the reconstructed 3D flow field of the Taylor droplet, we use the shear rates of the flow to visualize the location of high strain of the investigated Taylor droplets.

$$\begin{aligned} {\dot{\gamma }}^* =\frac{H_{\mathrm {ch}}}{8 u_{\mathrm {d}} } \left( \left( \left| \frac{\partial u}{\partial y}\right| +\left| \frac{\partial u}{\partial z}\right| \right) + \left( \left| \frac{\partial v}{\partial x}\right| +\left| \frac{\partial v}{\partial z}\right| \right) + \left( \left| \frac{\partial w}{\partial x}\right| +\left| \frac{\partial w}{\partial y}\right| \right) \right) \end{aligned}$$
(1)

As reference shear rate \({\dot{\gamma }}_{\mathrm {ref}}\) the wall shear rate of laminar flow in a circular channel with a diameter of \(H_{\mathrm {ch}}\) is applied: \({\dot{\gamma }}_{\mathrm {ref}} = \tfrac{8u_{\mathrm {0,tot}}}{H_{ch}} \approx \tfrac{8u_{\mathrm {d}}}{H_{\mathrm {ch}}}\). To enable a location of the strain maximum, we omit the directional influence of the derivatives and use absolute values only.

2.7 Vortex visualization applying the \(\lambda _2\)-criterion

We chose the \(\lambda _2\)-criterion (Jeong and Hussain 1995) to visualize the vortical structures in and around the Taylor droplet, which omits the effects from unsteady straining and viscosity. In their concept, a vortex is defined as a connected region with a negative second eigenvalue (\(\lambda _2 < 0\)), which corresponds to a vortex pressure minimum.

3 Results and discussion

In this section, we evaluate the µPIV measurement results of a 3D flow field relative to the interface of Taylor droplets in motion. At first, we detail the 3D velocity field on selected symmetry planes to introduce comparability to other studies and to visualize the hydrodynamic features. For this purpose, we refer to the vector component matrices u(x, y, z), v(x, y, z) and w(x, y, z) as U, V and W, respectively. Secondly, we analyze the shear rate distribution of the flow to indicate regions of maximum strain. Finally, the application of the \(\lambda _2\)-criterion allows visualization of the position of the wall-driven main vortices as well as the much less prominent compensating secondary vortex structures inside and outside the Taylor droplets.

If not stated otherwise, the velocity information is presented in a moving frame of reference relative to the droplet velocity \(u_{\mathrm {d}}\). The corresponding streamline representations of the flow are given either in symmetry planes or at locations, where the 2D reference section is not expected to exhibit velocities normal to the section (i.e., the droplet interface).

All representations of Taylor droplet interfaces in the following apply the analytical interface approximation of Mießner et al. (2019). Despite the large number of statistically incorporated Taylor droplets used in this study, we refer to a single Taylor droplet when describing and discussing the results.

3.1 2D2C velocity field at the center plane

At first, we discuss the normalized relative 2D2C field at the central symmetry plane of the Taylor droplet (Fig. 6) to assess the quality of the ensemble correlated multi-plane µPIV measurement. The droplets translate from left to right with a mean velocity of \(u_{\mathrm {d}} = {4.2}\, \hbox {mm}\, \hbox {s}^{-1}\). The droplet interface is indicated with the analytical shape approximation.

Fig. 6
figure 6

Exemplary ensemble correlated µPIV measurement of the normalized relative velocities inside and outside of a Taylor droplet at the x–y center plane of a square microchannel. The droplet translates from left to right. The yellow dot shows the downstream-pointing droplet tip at the origin coordinate system. The black contour line represents the interface position (Mießner et al. 2019). The velocity information is depicted in a moving frame of reference. a 2D2C vector plot of the measurement showing the flow field relative to the droplet movement (\(u_{\mathrm {d}} ={4.2}\, \hbox {mm}\, \hbox {s}^{-1}\)). For clarity, 3 out of 4 vectors in the x-direction are omitted. b Associated streamline plot of the relative flow field with the normalized stream-wise relative velocity being color-coded. c Associated streamline plot of the relative flow field with transverse velocities being color-coded. The color coding uses the radius \(r = \sqrt{x^2+y^2}\) as reference. Red represents flow away from the centerline (x-axis) in the positive radial direction, while the blue color indicates flow against the radial direction toward the centerline. The entire 3D2C data set is given in the supplementary material of this article, along with a plane-wise representation of the measured data

Full size image

Figure 6a depicts the vector field of the relative flow at the x–y symmetry plane. Only every fourth vector in the longitudinal direction is plotted to maintain clarity of the graph. The wall-driven main vortices are already recognizable in the vector plot.

In Fig. 6b, the corresponding streamline representation of the flow is shown in combination with a color coding of the normalized stream-wise relative velocity. The main vortices are evident inside and outside of the droplet. The maximum stream-wise absolute velocity at the center of the droplet’s main vortex is \(1.72 \, u_{\mathrm {d}}\). Also, the secondary vortex structures inside the front and the back cap of the droplet are captured clearly. The rear secondary vortex has roughly half the size of the secondary vortex in the front cap. The stagnation points at the cap tips are prominent.

The color coding in Fig. 6c shows the measured transverse velocity of the Taylor flow. We changed the reference axis of the flow direction to a radial representation \(r = \sqrt{x^2+y^2}\). The flow in the radial direction away from the x-axis is positive (red), while flow toward the center of the microchannel is negative (blue). The maximum radial velocities are a factor of about 3.5 times smaller than the droplet velocity \(u_{\mathrm {d}}\). Inside the secondary vortices at the front and back of the droplet, the magnitude of transverse velocities is a factor of 10 smaller than the droplet velocity \(u_{\mathrm {d}}\). Absolute values of the relative flow field inside the secondary vortices show that the local velocity magnitude decreases below \(0.005 \, u_{\mathrm {d}}\). Despite the large dynamic range of the measurement, the delicate flow structures of the secondary vortices are clearly present in the streamline plot.

At this point, we would like to refer the reader to the supplementary material, where the original 3D2C data set and a plane-wise representation of the measured data can be viewed.

3.2 3D3C velocity distribution

The streamlines on the Cartesian symmetry planes of the droplet allow the comparison of the quality of the measured velocity distribution (x–y center plane: U, V) to the reconstructed flow distribution (x–z center plane: U, W). The color-coded velocities of the flow relative to the interface are depicted in Fig. 7. The negative values indicate velocities directed against the chosen reference axis, while the positive values point into the positive direction. In addition, a streamline representation is given for each data slice based on the respective in-plane components. The droplets translate from left to right.

Fig. 7
figure 7

Measured (U, V) and reconstructed W velocity data of the normalized relative flow field in and around Taylor droplets. A streamline representation is given for each data slice based on the respective in-plane components. The droplets translate from left to right. The yellow dot shows the downstream pointing droplet tip at the origin coordinate system. The velocity information is depicted in a moving frame of reference. a The measured stream-wise relative velocity component is color-coded. The streamline representation of the flow field reveals symmetric features of the measurement in the x–y plane in comparison with the reconstruction in the x–z plane. b The radial transverse velocity component \(U_r(r,\theta ,x)\) is color-coded: blue pointing toward and red away from the x-axis (channel center). c Isosurfaces are added to b (\(\tfrac{U_r}{u_{\mathrm {d}}} = 0.175\)) to visualize the transverse velocities

Full size image

The stream-wise velocity components are emphasized in Fig. 7a; thus, solely measured information is color-coded. Despite the difference in spatial resolution between data in the x–y plane (\(\varDelta x = \varDelta y ={1.8}\, \upmu \hbox {m}\)) and the resolution in the z-direction (\(\varDelta z = {5.4}\, \upmu \hbox {m}\)), the symmetry of the measurement with respect to the x-axis is obvious. Approaching the droplet rear cap from the tailing slug, the fluid at the channel centerline decelerates in the proximity of the stagnation point. Through the entire rear secondary vortex, the liquid moves slowly with almost the x-velocity of the interface \(u_{\mathrm {d}}\). Further, along the x-axis the fluid accelerates, when it enters the wall-driven main vortex and decelerates toward the frontal secondary structure. Beyond the frontal stagnation point the outer fluid accelerates again into the wall-driven main recirculation. In the proximity of the wall, proceeding from the droplet front to its back cap, the fluid decelerates at the on- and off-set of the wall film.

Based on the reconstructed third velocity component W in the z-direction (explained in Sect. 2.4) and the measured transverse velocity V, we show the normalized radial transverse velocity \(U_r(r,\theta ,x)\) in Fig. 7b. The color coding labels velocities blue when pointing toward and red when pointing away from the x-axis (channel center). The reconstruction of \(W=U_r\) in the x–z plane with \(y=0\), where \(V(y=0) = 0\), is in good agreement with the measured transverse velocities \(V=U_r\) in the x–y plane with \(z=0\), where \(W(z=0)=0\). Also, the amplitude of both positive and negative transverse velocities inside the droplet is similar. However, the velocity development from the stagnation points toward the upper channel wall shows that the reconstruction overestimates the z-component.

In general, the flow features along the z-axis are less frail in comparison with the x–y symmetry plane due to the lower spatial resolution. Additionally, the streamline representation of the flow field reveals the highly symmetric features of the measurement in the x–y plane in comparison with the reconstruction in the x–z plane.

Isosurfaces of the transverse velocities are added (\(\tfrac{U_r}{u_{\mathrm {d}}} = 0.175\)) to visualize the symmetry of the 3D spatial distribution of the transverse velocities (Fig. 7c). The features of the isosurfaces are generally symmetrical. However, their relative position at the channel top compared to the location with respect to the sidewalls reveals slight asymmetry. The amplitude of the reconstructed velocities in the z-direction is too pronounced: While the isosurfaces of the disperse phase remain inside the droplet in the proximity of the sidewall, they clearly breach the interface toward the top wall. This corresponds to the considerations on the volume flow described in Sect. 2.5.

3.3 Relative velocity field on the interface of the Taylor droplet

The fluid motion that is directed perpendicular to a droplet interface causes the interface to deviate from its shape at rest (Mießner et al. 2019). A clean interface cannot resist shear forces; thus, tangential velocity components cause the fluid at the interface to move: Momentum coupling across the interface transfers kinetic energy into the droplet phase. Thus, the flow of the continuous phase drives the velocity field inside the droplet.

Figure 8a shows the vector field of the interpolated velocity distribution on the interface. The top wall film is covered with a white area to mark a lack of available data. The velocity magnitude on the front and back caps is smaller than the interface movement in the proximity of the wall films. The interface movement in the gutter region is directed from the front to the back of the droplet. The velocities along the gutter interface range between the fast wall movement and the slow movements on the caps.

Fig. 8
figure 8

Interpolation of the 3D velocity field on the interface of the Taylor droplet. The velocity information is depicted in a moving frame of reference. a In-plane 2D vector representation of the velocity field on the mobile droplet interface. The sampling is reduced in the x–y-direction to improve the clarity of the graph. b Streamline plot of the in-plane velocity distribution on the Taylor droplet interface. The color coding indicates the normalized stream-wise relative velocity of the flow. In the droplet frame of motion, the walls move in the negative x-direction. At the droplet top, no data are available. Thus, the section is covered with a corresponding white surface. c In addition to the streamlines on the droplet interface, the color coding shows the radius-dependent transverse velocity component (red—away from the x-axis; blue—toward the x-axis). d Rear view d and front view e of the droplet indicating the same quantities as in c

Full size image

In Fig. 8b, the direction of the interface motion is given as a streamline representation. The color coding indicates the normalized relative velocity component in the stream-wise direction. Fluid moving in the positive x-direction relative to the velocity of the droplet is shown in red, while blue indicates motion against the downstream direction. White color shows interface regions that do not move with respect to the overall translation of the droplet. As expected, two stagnation points are identified: One is located at the front tip and the other at the back tip of the droplet. In addition, two annular stagnation regions can be seen, where the forward moving core of the droplet flow separates from the backward-directed motion of the wall-driven interface.

The radius-dependent transverse velocities of the interface movement are shown in Fig. 8c as color coding. Red shows the fluid motion outward to the wall, while blue indicates the inward movement toward the x-axis. In combination with the streamlines, the stagnation points at the droplet front and back tips can be seen clearly. At the front tip, fluid elements of the interface converge into the stagnation point. At the back tip, the diverging stagnation point shows that a new interface is generated.

The annular stagnation region at the on- and offset of the wall films (and the gutters) is clearly recognizable in the transverse component of the interface motion. At the droplet front, fluid elements below a certain distance from the channel center are pulled inward to the stagnation point. Fluid elements beyond that distance are directed outward and are either transported into the wall films or through the gutters. At the end of the gutter and of the wall films, the interfacial fluid is transported inward toward the rear annular stagnation region. Here, it meets with the outward-directed interface motion, which is driven by the rear main vortex.

The annular stagnation regions can be observed best in a rear view (Fig. 8d) and a front view (Fig. 8e) of the Taylor droplet. The color coding shows the transverse velocity component on the interface. Figure 8d, e does not show the perfect symmetry which we would expect. We attribute the deviation of the transverse velocity in the z-direction to the lower measurement resolution and to the reconstruction accuracy of the out-of-plane velocity component. The diverging stagnation point of the rear cap tip (Fig. 8d) is complemented with a converging annular stagnation region. We find a small vortex structure where the fluid elements meet, that leave the wall film and the gutter. Symmetry suggests that a similar structure should be present at the corresponding location on the wall in the y-direction. However, the spatial resolution in the z-direction does not permit the detection of any small filaments. In the front view (Fig. 8e) of the droplet interface motion, a divergent annular stagnation region produces a new interface, which converges into the frontal stagnation point.

A comparison of the relative radial location of the two annular stagnation regions shows that both of them follow the shape of the channel cross section. However, the frontal stagnation region reaches in an outward direction further into the gutter region and results in a pillow shape. In contrast, the rear stagnation region is shifted inward toward the center axis of the microchannel and maintains a more oval shape.

3.4 Stagnation regions

The 3D distribution of stagnation regions inside and outside of the Taylor droplet (Fig. 9) is investigated using the relative flow field. We compute its absolute magnitude, normalize the scalar values with the droplet velocity (\(\sqrt{( U - u_{\mathrm {d}} )^2 + V^2 + W^2} /u_{\mathrm {d}}\)) and project the result on the interface and the x–y and x–z symmetry planes. In addition, isosurfaces of the relative velocity magnitude at a dimensionless level of 0.075 indicate stagnation regions. The threshold corresponds roughly to a fraction of less than 1% of the present continuous phase kinetic energy. It is chosen such that all stagnation regions are clearly visible.

Fig. 9
figure 9

Visualization of stagnation regions using isosurfaces (Threshold: 0.075). The depicted scalar quantity is the absolute value of the relative velocity field. The velocity field information is depicted in a moving frame of reference. Rear view (a), front view (b) and elevated side view (c) visualize the 3D location of the stagnation regions with respect to the droplet interface. In a, b and c the streamlines on the interface indicate the motion of the droplet surface. d reveals according stagnation regions inside the droplet and the location of the core of the main disperse phase vortex. The interface shape at the droplet top film is covered in white, since no data is available that can be projected to this location

Full size image

At the back (Fig. 9a) and at the front (Fig. 9b) of the droplet, the stagnation region breaches the interface and indicates its location with respect to the droplet shape. At the droplet back the region is connected, while at the front, the stagnation point at the droplet tip is separated from an annular stagnation region. The outer shape of both stagnation regions clearly follow the cross section shape of the microchannel.

To observe the relative location of the stagnation regions with respect to each other, Fig. 9c shows an elevated side perspective. At the location of the annular stagnation region, the outer vortices attach to (droplet front) or detach from (droplet back) the interface. Here, the hydrodynamic forces of the flow interact strongly with the interface to deform the droplet (Helmers et al. 2019a). The corner of both outer stagnation regions marks the entrance of the gutter, where continuous liquid bypasses the droplet. The streamline pattern of the interface motion coincides with the visualization of the absolute magnitude of the relative velocity.

Without the interface, one can observe the entire extent of the stagnation regions (Fig. 9d). The isosurfaces cover a major part of the secondary vortices inside the droplet caps. Additionally, the annular core structure of the wall-driven main vortex of the disperse is prominent, since it is a region of low relative velocity compared to the motion of the Taylor droplet. The position agreement of the main vortex center on the x–z symmetry plane and the location of 3D stagnation region at the vortex core demonstrates the symmetry of the 3D flow field.

3.5 Shear rate distribution

The shear rate is related to the magnitude of the strain. We use its normalized absolute amplitude to establish a scalar quantity to visualize the strain distribution in and around the Taylor droplet as well as to locate the strain maxima (Fig. 10).

Fig. 10
figure 10

Absolute values of the normalized shear rate inside and around Taylor droplet indicating strain. The velocity information is depicted in a moving frame of reference. a Normalized shear rate distribution on the symmetry planes (x–y and x–z directions) of the flow. b Isosurfaces (dimensionless threshold level: 1.5) to expose the location of maximum strain. The maximum strain is located at gutter and wall film entrances and respective exits. The interface model is artificially closed at the top film

Full size image

The shear rate distribution on the symmetry planes in the x–y and x–z directions show the expected general increase from the microchannel center toward the walls (Fig. 10a). The wall shear rate decreases in the proximity of the cap interface, since at these stagnation regions the velocities and their gradients are low. At the longitudinal center of the main vortices, the relative wall motion that drives the flow causes a clear shear rate increase. A comparison of the shear rate on the two symmetry planes reveals a general similarity. However, the errors introduced by the reconstruction (above the cap stagnation points in the z-direction) and the noise increase due to the derivatives have a recognizable impact. In addition, delicate structures on the x–y plane cannot be observed in the perpendicular plane.

The normalized shear rate information at the interface position and isosurfaces of an elevated shear rate (dimensionless threshold level: 1.5) enables analysis of the impact of the strain location on the Taylor droplet shape (Fig. 10b). At the droplet front cap, the strain location coincides with the position of pressure buildup. Fluid decelerates and the continuous phase either evades the interface and converges into the gutter and the wall films or the flow is forced to reverse with respect to the droplet translation. The position of the shear rate maximum overlaps with the corner of the annular stagnation region, which can be identified in Fig. 9c. The pressure increase at the stagnation region squeezes and elongates the front cap in the flow direction from its spherical shape at rest.

At the droplet back cap, the strain location coincides with the position of pressure release. The continuous phase accelerates and diverges from the wall films and the gutters until it merges with the reversing flow of the rear main vortex. The location of the shear rate maximum at the back coincides with the corner position of the annular stagnation region (Fig. 9c). The interface is being pulled toward the wall. In combination with the pressure buildup at the rear stagnation point, a compressed deformation in the main flow direction is achieved with respect to the undeformed static droplet shape.

The observation of maximum shear in the continuous phase at the entrances and exits of the gutters leads to the recommendation to place shear-sensitive reactants or products in the dispersed phase.

3.6 Vortex strength: \(\lambda _2\)-criterion

The \(\lambda _2\)-criterion corresponds to the pressure minimum on a planar section through the flow and is applied to represent vortex core geometry (Jeong and Hussain 1995). A vortex is defined as a connected region, where \(\lambda _{2}\) is negative. Hence, only negative values are depicted in Fig. 11.

Fig. 11
figure 11

Application of the \(\lambda _2\)-criterion to visualize the vortices inside and outside of the Taylor droplet. Solely negative values are depicted, since only connected negative regions represent the pressure minima of vortices. The presented velocity information is given  in a moving frame of reference. a Scalar representation of the \(\lambda _2\)-criterion on the symmetry planes of the Taylor flow (x–y and x–z plane through the origin): The derivatives of the velocity reconstruction introduce a significant amount of noise. b The fluid motion on the interface causes the inner circulation of the Taylor droplet. The streamline pattern of the interface motion is depicted separately to improve the clarity of Fig. 11c, d. c Isosurface visualization of the vortices inside the disperse phase only: The main vortex is about four times stronger compared to the front vortex, which in turn is about three times stronger than the rear secondary structure. Derivative-induced noise is clearly visible at the main vortex in the proximity of the top wall. The threshold of the main vortex (blue isosurface) is set to − 1.28 to follow the outline of the vortex in the x–y symmetry plane. For the secondary vortices, we chose different threshold levels for the front (− 0.24, magenta isosurface) and the back of the droplet (− 0.08, green isosurface) to visualize their 3D extent inside their respective caps. d\(\lambda _2\) representation of the continuous phase flow is added: The recognition of coherent structures becomes increasingly difficult due to the noisy data at the wall

Full size image

Observation of two sections through the scalar field of the \(\lambda _2\)-criterion on the symmetry planes of the Taylor flow (x–y and x–z directions) (Fig. 11a) enable us to detect the main features and to judge the data quality. On the x–z plane, the main wall-driven vortices are as clearly recognizable as they are on the delicately featured x–y plane. However, the velocity derivatives in the z-direction amplify the noise of the reconstructed velocity z-component.

The fluid motion on the interface causes the inner circulation of the Taylor droplet. The streamline pattern of the interface motion is depicted separately from the vortex visualization to maintain clarity (Fig. 11b). We use isosurfaces only inside the disperse phase to visualize both: The large main vortex and the smaller counter rotating secondary structures of the front and back of the droplet caps can be clearly distinguished (Fig. 11c). The main vortex (blue isosurface) is not strictly toroidal as one might expect from Taylor flows in circular tubes, but its structure reflects the square channel geometry. The extent of the main vortex in the plane of the cross section diagonal and along the main flow axis is decreased toward the droplet front. Nonetheless, the main recirculation maintains its symmetry with respect to the channel dimensions. The noise, attributed to the squared and summed derivatives of the reconstructed velocity components, is noticeable in the proximity of the top wall of the channel.

The location of the secondary vortices (green and magenta isosurface) that compensate for the rotation direction of the main recirculation coincide with the stagnation region identified in Fig. 9. Despite their low values of \(\lambda _2\), both structures are clearly detectable and carry a channel-related symmetry. In addition, the shape and orientation of the secondary vortices match the location of the streamline pattern on the droplet interface (Fig. 11b).

In the droplet frame of motion, the wall movement directly drives the main vortices in the slug between the droplets. It also moves the separating film between the wall and the droplet. The flow inside the droplet is generated by the momentum transfer across the interface, and therefore the velocity on the interface drives both the main and the secondary vortices inside the droplet. The magnitudes of the interface velocities determine the spatial extension of vortices. While the main vortex is powered by the strong stream-wise flow of the wall film, the secondary vortices are driven by transverse velocities, which are an order of magnitude smaller.

Adding the \(\lambda _2\) vortex detection of the continuous phase to the graph, the presence of the main vortex of the adjacent slug can only be anticipated (Fig. 11d), yellow isosurface). Despite choosing the identical threshold as for the main vortex in the droplet (− 1.28), the noise of the reconstructed data set prevents a clear detection of its structure.

4 Conclusion

We conducted a µPIV study on the 3D2C velocity distribution of Taylor droplets moving in a square horizontal microchannel. We reconstructed the third velocity component and verified the result based on the volume flow in stream-wise and transverse flow directions. The Taylor droplets dilate the PDMS channel walls and cause the stream-wise volume flow to be underestimated up to \(-4.7\%\). The mean deviation of reconstructed velocity in the z-direction is found to be \(\pm 2.7\%\), while the maximum deviation amounts to − 6% of \(Q_{\mathrm {tot}}\). In the y-direction, the deviation of the ensemble correlation µPIV is quantified to reside below \(\pm 0.1 \%\).

The completed 3D3C velocity field is used to investigate the stagnation regions of the 3D flow. In combination with the interface approximation of Mießner et al. (2019), the location of the stagnation regions reveals expectable stagnation points at the cap tips of the Taylor droplet. In addition, we visualized an annular stagnation region closer to the walls, where the main outer vortices attach to or detach from the interface. It is at this location that the hydrodynamic forces applied on the interface cause the droplet deformation.

We located the maximum shear rate at the entrances and exits of the wall films and the corner flow (gutter) bypassing the Taylor droplet. This straining location macroscopically corresponds to the cap position of the Taylor droplets.

For the first time, an experimental data set allows to visualize the streamlines of the velocity distribution on the interface of a Taylor droplet and to directly relate it to the main and secondary vortices of the velocity field. We found a small vortex structure in the stagnation region at the back of the droplet, where the fluid elements that leave the wall film and the gutter meet.

The evaluation methods which have been developed for this study allow to investigate laminar quasi-stationary and pulsating flows. Determining the shear rate distribution helps to effectively handle fragile educts and/or products in single phase as well as in multiphase microfluidic flows, e.g., organ-on-a-chip applications and cell handling applications using droplets (Cui and Wang 2019).

5 Outlook

To further exploit the 3D measurement, we will proceed with deduced quantities of higher order. In a follow-up study, we plan to use this data set in the Navier–Stokes equation to investigate the flow-induced pressure contribution to the deformation of the Taylor droplet from its static shape at rest. We aim to make use of the interface curvature information delivered by the approximate interface model (Mießner et al. 2019) to derive Laplace pressure and geometry-based deformation–pressure distribution. After a qualitative comparison of both results, we will propose two methods to estimate the pressure drop of the Taylor droplet and to estimate the driving pressure gradient for the bypass flow through the gutter of the Taylor droplet. The latter gives rise to the velocity difference between the mean flow and the Taylor droplet (Helmers et al. 2019b).