Modeling the excess velocity of low-viscous Taylor droplets in square microchannels

TL;DR

This paper presents a new model to predict the excess velocity of Taylor droplets in square microchannels, linking droplet shape deformation, bypass flow, and flow parameters, validated with experimental data.

Contribution

The work introduces a novel model for droplet excess velocity considering cap deformation and bypass flow, validated with empirical data and optimized via metaheuristic methods.

Findings

Model accurately predicts droplet excess velocity.

Parameters like bypass length and viscosity ratio are key.

Validated against published experimental data.

Abstract

Microscopic multiphase flows have gained broad interest due to their capability to transfer processes into new operational windows and achieving significant process intensification. However, the hydrodynamic behavior of Taylor droplets is not yet entirely understood. In this work, we introduce a model to determine the excess velocity of Taylor droplets in square microchannels. This velocity difference between the droplet and the total superficial velocity of the flow has a direct influence on the droplet residence time and is linked to the pressure drop. Since the droplet does not occupy the entire channel cross section, it enables the continuous phase to bypass the droplet through the corners. A consideration of the continuity equation generally relates the excess velocity to the mean flow velocity. We base the quantification of the bypass flow on a correlation for the droplet cap deformation from its static shape. The cap deformation reveals the forces of the flowing liquids exerted onto the interface and allows estimating the local driving pressure gradient for the bypass flow. The characterizing parameters are identified as the bypass length, the wall film thickness, the viscosity ratio between both phases and the Ca-number. The proposed model is adapted with a stochastic, metaheuristic optimization approach based on high-speed camera measurements. In addition, our model is successfully verified with published empirical data.

Tables

Table 1. Table 1: Boundaries of input values for optimization

parameter	lower boundary	upper boundary
$m_{cap,f}$	1.00	9.90
$c_{cap,f}$	0.40	1.50
$n_{cap,f}$	1.00	1.005
$m_{cap,b}$	-2.50	-1.00
$c_{cap,b}$	0.30	0.75
$n_{cap,b}$	0.995	1.00
$m_{\beta }$	0.00	10.00
$c_{\beta }$	0.50	1.5
$n_{\beta }$	0.50	20

Table 2. Table 2: Weight factors ω i subscript 𝜔 𝑖 \omega_{i} of the loss-function for GA optimization

${\omega }_1$	${\omega }_2$	${\omega }_3$
3.0	5.0	3.7

Table 3. Table 3: Properties of used solver algorithms

Genetic algorithm
Population size	200 individuals
Creation function	Random feasable population
Scaling function	Ranking
Selection function	Stochastic uniform
Mutation function	Adaptive feasible
Crossover function	Scattered
Pattern Search algorithm
Search method	Latin Hypercube
Poll method	Complete poll

Table 4. Table 4: Values for the fit-functions for droplet shape k i subscript 𝑘 𝑖 k_{i} and resistance β 𝛽 \beta

Target value	$m_i$	$c_i$	$n_i$
$k_f$	4.8761	0.7465	1.0021
$k_b$	-1.6967	0.4745	0.9980
$\beta$	6.1280	1.2105	1.4541

Equations

$$Ca=\frac{u_0{\eta }_c}{\sigma }$$

$$u_0=\frac{Q_0}{A_{ch}}=\frac{Q_d+Q_c}{W\cdot H}=\frac{Q_d+Q_c}{H^2\cdot ar}$$

$$Re=\frac{u_{0}d{\rho }_c}{{\eta }_c}$$

$$\lambda =\frac{{\eta }_d}{{\eta }_c}$$

$$Oh=\frac{{\eta }_c}{\sqrt{d{\rho }_c\sigma }}=\sqrt{\frac{Ca}{Re}}$$

$$u_{slip}=\frac{u_d-u_0}{u_d}=1-\frac{u_0}{u_d}$$

$$u_{ex}=\frac{u_d-u_0}{u_0}=\frac{u_d}{u_0}-1$$

$$Q_0=Q_d+Q_c=Q_d+\left(Q_s+Q_g+Q_f\right)$$

$$Q_{rel}=u_{rel}\cdot A_{ch}$$

$$Q_0-Q_{rel}=Q_d+Q_g+Q_f-Q_{rel}$$

$$A_g$$

$$A_f$$

$$A_d$$

$$u_0{A}_{ch}-u_{rel}{A}_{ch}=\overline{u_d}{A}_d+\overline{u_g}{A}_g+\overline{u_f}{A}_f-u_{rel}{A}_g-u_{rel}{A}_d-u_{rel}{A}_f$$

$$\overline{u_d}=u_d\stackrel{!}{=}{u}_{rel}$$

$$\overline{u_f}\approx 0$$

$$(u_0-u_d)A_{ch}=\overline{u_g^{rel}}{A}_g+(0-u_d)A_f$$

$$\frac{u_d}{u_0}-1=u_{ex}=-\frac{\overline{u_g^{rel}}{A}_g}{u_0{A}_{ch}}+(u_{ex}+1)\frac{A_f}{A_{ch}}$$

$$u_{ex}=-\frac{Q_g^{rel}}{Q_0}+\frac{A_f}{A_{ch}}$$

$$k_{g,i}=\frac{R_{g,i}}{h}=\frac{R_{g,i}}{H-2\delta }$$

$$\overline{R_g}$$

$$=\frac{(H-2\delta )(k_{g,f}+k_{g,b})}{2}$$

$$k_{c,f}=\frac{a_f}{b_f}$$

$$k_{c,b}=\frac{a_b}{b_b}$$

$$k_{c,i}=m_{cap}\cdot C{a}^{c_{cap}}+n_{cap}$$

$$\Delta p_{g,f}=\sigma (\frac{1}{R_{g,f}^{stat}}-\frac{1}{R_{g,f}})$$

$$\Delta p_{g,b}=\sigma (\frac{1}{R_{g,b}^{stat}}-\frac{1}{R_{g,b}})$$

$$\Delta p_{LP,fb}$$

$$=\sigma (\frac{1}{R_{g,f}}-\frac{1}{R_{g,b}})$$

$$\Delta p_{LP,fb}$$

$$=(\frac{\sigma }{H})\frac{1}{(1-\frac{2\delta }{H})}(\frac{k_{g,b}-k_{g,f}}{k_{g,f}\cdot k_{g,b}})$$

$$Q_g^{rel}=\frac{1}{\Omega }\Delta p_{LP,fb}$$

Full text

\pubvolume

xx \issuenumtba \articlenumbertba

\historyReceived: tba; Accepted: tba; Published: tba

\TitleModeling the excess velocity of low-viscous Taylor droplets in square microchannels \AuthorThorben Helmers1\orcidA*, Philip Kemper 2, Jorg Thöming 2\orcidC and Ulrich Mießner 1\orcidD \AuthorNamesThorben Helmers, Philip Kemper, Jorg Thöming and Ulrich Mießner

\corresCorrespondence: helmers@uvt.uni-bremen.de; Tel.: +49-421-218-63337

1 Introduction

Microscopic multiphase flows facilitate a wide field of possible applications since they provide short diffusion layers within the flow structures. This enables high mass and heat transfer rates (Haase et al., 2016; Sattari-Najafabadi et al., 2018) for several applications ranging from extraction (Kralj et al., 2007) and multiphase catalyst reactions (Kobayashi et al., 2006) to improved unit operations like mixing tasks (Wong et al., 2004). The distinct features allow performing reactions at new process windows with fewer hazards or higher selectivity (Sun et al., 2016). The specific flow conditions can furthermore serve for cell isolation (Chen et al., 2014), genetic analysis (Hosokawa et al., 2017) and reaction screening in a droplet chain (Jin et al., 2018). In contrast to large scale multiphase flows, microscopic flows are much easier to predict as there are no complex interactions such as swarm turbulence (Kück et al., 2018; Mießner et al., 2017) commonly found in bubble columns. In fact, the reproducibility of e.g. Taylor flows is a key for the application of microscopic multiphase flows (Kreutzer et al., 2005).

In the Taylor flow regime, the disperse phase is separated from the wall by a thin wall film and does not fill the entire cross-section of the microchannel. The remaining space between the droplet and the microchannel corners is occupied by the continuous phase, which are referred to as gutters (Vansteijn et al., 2008). The droplets are typically longer than the channel diameter, which leads to separated elongated disperse phase instances. The continuous phase segments between the droplets are called slugs. Taylor flows are mainly established in circular capillaries or rectangular microchannels whereas especially the hydrodynamics of moving Taylor droplets in circular capillaries and the role of the thin wall film have been intensively studied (Abiev et al., 2017).

Chemical reactions on the microscale are often performed in monolith reactors as a catalyst support. Within these reactors, a high number of parallelized channels with hexagonal or square channel cross-section offer a high specific reaction area for wall placed catalysts at small wall thickness. This results in a better heat transfer through the walls and better mechanical stability than circular capillaries (Gascon et al., 2015). For process control and stabilization, as well as precise reactor design, knowledge of the underlying fluidic terms is crucial. The high grade of parallelity complicates prediction of the hydrodynamics and the resulting pressure drop (Kreutzer et al., 2005; Huerre et al., 2014).

Besides disperse phase size distribution and formation frequency (Fu and Ma, 2015), the actual droplet velocity is essential for the droplet residence time in the reactor. It determines the contact time of the educts and influences the pressure drop of the reactor (Ładosz and Rudolf von Rohr, 2018). Within a parallel reactor, an exact knowledge of the pressure drop is especially necessary since a steady educt supply for each individual single reactor is needed to ensure stable and efficient working conditions (Schubert et al., 2016).

Several publications have dealt with the droplet velocity in rectangular capillaries and observed a droplet velocity mostly faster than the superficial velocity (see Sec. 3). For flows within circular capillaries, where only a thin wall film is present, this velocity difference is well understood (Abiev, 2017), while for rectangular microchannel a variety of explanations exist, that mostly correlate the relations from measurements (Jose and Cubaud, 2014). This complicates the transfer of results to other flow applications or altered process parameters since local and instantaneous hydrodynamic parameters are mostly not taken into account by the models and correlations.

This work aims to establish a model to determine the droplet velocity from the actual flow conditions: e.g. droplet length, material properties, and the $Ca$-number. In a first step, we develop a concept for the relative droplet velocity, which bases this velocity on extrinsic parameters, allowing symmetrical scaling. From this concept, we identify the bypass flow through the gutters as well as the film-thickness as the prominent parameters for the excess velocity.

In the next step, we develop a model that uses the local surface curvature of the gutters to retrieve the local pressure at the entrance and outlet of the gutter’s corner as a driving force. The bypassing gutter flow is calculated based on the counterplay of this driving force, the gutter length and a viscosity correlated resistance factor $\beta$. The local droplet curvature at the gutter entrances is derived with an analytical interface shape approximation (Mießner et al., 2019) and a correlation for the droplet cap curvature based on the $Ca$-number from our previous work (Helmers et al., 2019). The model is successfully validated by high-speed camera measurements.

2 Hydrodynamic fundamentals of Taylor Flows

The Taylor flow regime in rectangular micro-channels is mainly influenced by surface tension forces rather than inertia forces. In the Taylor flow regime, a droplet fills nearly the whole cross-section of a hydrodynamic channel, while the continuous phase occupies the gutters and a thin wall film. The droplets are divided by the slugs, consisting only of liquid from the continuous phase.

The interaction between interfacial and viscous forces is described by the Capillary number $Ca$.

[TABLE]

Herein, $\eta_{c}$ represents the dynamic viscosity of the continuous phase, $\sigma$ the interfacial tension between both phases and $u_{0}$ the total superficial flow velocity:

[TABLE]

The superficial velocity $u_{0}$ is based on the volume flow of the disperse ($Q_{d}$) and continuous ($Q_{c}$) phase as well as the microchannel’s cross-sectional area $A_{ch}$ that calculates from the channel width $W$ and channel height $H$ or respectively the channel aspect ratio $ar=\frac{W}{H}$.

For energetic considerations, knowledge of the ratio between inertia and viscous forces is of importance. The Reynolds-number

[TABLE]

is used to describe the ratio between inertia and viscous forces of the continuous phase with $\rho_{c}$ being the continuous phase density and $\eta_{c}$ the continuous phase viscosity. Additionally, the viscosity ratio

[TABLE]

between both phases has a significant influence on the pressure drop and the velocity of a droplet (Direito et al., 2017; Shams et al., 2018; Ładosz and Rudolf von Rohr, 2018) since it indicates the momentum coupling into the disperse phase. Please note, that the definition of $\lambda$ differs within the mentioned publications.

Considering the overall classification of the applied flow system, the material properties of both fluid phases are of importance. The Ohnesorge-number $Oh$ describes the most prominent material properties for droplets being formed or dispersed (Ohnesorge, 1936). It characterizes the fluidic system independent of current flow or forces and is mostly used when working with surfactants to manipulate the flow properties. Here, we use the $Oh$-number to characterize the continuous phase.

[TABLE]

In many applications, Taylor droplets in rectangular microchannels move with a velocity different from the superficial velocity, because the droplet does not fill the entire channel cross section and continuous phase can bypass the droplet through the gutters. An early description is given by Wong et al. (1995), who dewscribes these phenomena analytically and declare two possible regimes.

In the first regime, the fluid in the gutters moves slower than the droplet, dissipating kinetic energy. For the gutters, Abiev (2017) reported an inversed pressure gradient, indicated by the local surface curvature. For circular microchannels, this results in an inverse flow of the wall film and the droplet moves faster than the superficial velocity.

In the second flow regime, which holds true for long and highly viscous droplets (Wong et al., 1995; Shams et al., 2018; Rao and Wong, 2018), more energy is dissipated through the larger wall film area and through viscous dissipation within the droplets. Consequently, the flow in the gutters moves from the droplet back to the front. Thus, the droplet moves slower than the superficial velocity.

For both regimes, the thin wall film resists the motion of the droplet, resulting in a difference in pressure with higher value at the back and a lower value at the front of the droplet. The transition between both regimes is described by a critical flow rate and depends on the droplet length and the channel aspect ratio (Wong et al., 1995).

Jakiela et al. (2011) focus experimentally on the influence of the momentum coupling between both phases represented by $\lambda$ and also reveal a dependence of the droplet velocity on the droplet length. For short low viscous droplets ($\lambda<1$ ), the droplets move faster than the superficial velocity and the droplet behavior is assigned to the first flow regime. Highly viscous droplets ($\lambda>1$ ) move either faster or slower than the superficial velocity, depending strongly on the droplet length.

3 Concept of Excess Velocity

Based on the instantaneous droplet velocity $u_{d}$, that is evident and directly measurable via optical or electrical measurement techniques (Kalantarifard et al., 2018), different explanatory approaches for the deviation of superficial velocity $u_{0}$ and droplet velocity $u_{d}$ have been reported.

Liu et al. (2005) define this relative difference as slipping velocity $u_{slip}$, Howard and Walsh (2013) as relative drift velocity $u_{drift}$, Angeli and Gavriilidis (2008) as relative bubble velocity $u_{rel}$ and Abadie et al. (2012) as dimensionless droplet velocity. Jakiela et al. (2011) focus directly on the ratio of $m=\frac{u_{d}}{u_{0}}$ and name this quotient droplet mobility according to Bretherton (1961). For this work, we like to summarize these approaches as a slipping velocity:

[TABLE]

In those concepts, the desired quantity is scaled with an intrinsic value such as the instantaneous droplet velocity, which leads to normalization effects as values $u_{slip}<1$ and $u_{slip}>1$ are not normalized symmetrically (Fig. 1). This behavior has to be taken into account when experimental or simulative data is interpreted.

In our approach, we scale the velocity difference with the superficial velocity as an extrinsic property and, staying in the term of extrinsic denomination, define it as an excess velocity $u_{ex}$

[TABLE]

In this manner, $u_{ex}$ results values around [math] for droplet velocities equal the superficial velocity (plug flow), while positive and negative values indicate droplets, which are moving respectively faster or slower than the superficial velocity.

The advantage of this extrinsic concept stands out in a comparison of both approaches (Fig 1). The first shown intrinsic concept (slip velocity) leads to an asymmetrical scaling behavior especially for droplets with $u_{d}<u_{0}$, since $u_{slip}<0$ decreases stronger than $u_{slip}>0$ would increase. For applications like balancing or process modeling, a linear behavior is mandatory, to prevent an additional bias of the modeled quantities towards any direction $u_{d}<u_{0}\lor u_{d}>u_{0}$.

A description for the excess velocity can be derived from the volume flows around moving droplets (Eq. 8). The continuity equation describes the interrelation between the gutter flow and the outer driving flows and it delivers the relation between the total flow ($Q_{0}$) and the volume flow fractions of the disperse ($Q_{d}$) as well as the continuous phase ($Q_{c}$). Considering the unit cell of a single slug and an adjoining droplet, the continuous phase parts into the volume flow of the slug ($Q_{s}$), through the gutters ($Q_{g}$) and through the wall film ($Q_{f}$):

[TABLE]

If we depict the flows at a moving Taylor droplet (Fig. 2 a) ) and introduce a stationary control surface $\Gamma$ (Fig. 2 b) ), velocities can be retrieved from the balances, while two differing flow states are possible: In the case of a droplet passing $\Gamma$, the slug volume flow through the control surface is $Q_{s}=0$ (since there is no slug present). For a slug passing $\Gamma$, only the slug volume $Q_{s}$ is present. As one can see, the use of a stationary point of view leads to instationary terms within the balances.

If the coordinate system is changed from a Lagrangian to an Eulerian system by moving the control surface with an arbitrary velocity $u_{rel}$, the balances become stationary and relative velocities become visible (Fig. 2 c) ). For this moving coordinate system, an additional volume flow $Q_{rel}$ adds to the balances, that results from the transformation of the coordinate system (Fig. 2 d):

[TABLE]

This changes Eq. 8 to:

[TABLE]

With knowledge of the specific areas for each distinct volume flow rate, the averaged velocities can be calculated. Following 3 a) + c) we conclude for the gutter area $A_{g}$ of all four gutters

[TABLE]

with $\delta$ denoting the wall film thickness and $\overline{R_{g}}$ representing the mean dynamic gutter radius, which is described later in Sec. 4. For the cross-sectional film area $A_{f}$ we state with the aspect ratio $ar=WH^{-1}$

[TABLE]

The droplet cross-sectional area $A_{d}$ is delivered combining ${A_{g}}$ and $A_{f}$:

[TABLE]

With the cross-sectional areas of the droplet, the gutter and the wall film, the volume flow rates can be rearranged to area-averaged velocities and Eq. 10 becomes

[TABLE]

Herein $\overline{u_{d}}$, $\overline{u_{g}}$ and $\overline{u_{f}}$ are area-averaged velocities of the droplet, gutter and film.

For the droplet ($Q_{d}=\overline{u_{d}}A_{d}-u_{rel}A_{d}$) and the film volume flow ($Q_{f}=\overline{u_{f}}A_{f}-u_{rel}A_{f}$), the transition velocity $u_{rel}$ equals the stagnation point velocity $u_{d}$, since we assume incompressibility, mass conservation, a stagnant film and a stationary droplet shape (Luis A. M. Rocha and João M. Miranda and Joao B. L. M. Campos, 2017)

[TABLE]

Additionally we assume the averaged velocity in the thin wall films to be insignificant for $Ca<0.2$

[TABLE]

Therefore Eq. 14 simplifies to

[TABLE]

Herein $\overline{u_{g}^{rel}}=\overline{u_{g}}-u_{d}$ represents the relative gutter velocity, which dissipates flow energy within the gutter. Simplification and combination with Eq. 7 leads to

[TABLE]

We neglect the terms, that are small of higher order (see Appendix A) and retrieve an expression for the excess velocity:

[TABLE]

The relative gutter volume flow ($Q_{g}^{rel}$) and the cross-sectional area of the wall film ($A_{f}$) are the most prominent influencing quantities for the excess velocity. Thus our proposed model aims to especially determine these quantities.

4 Model Specification

The considerations of the previous section identify the volume flows significantly determining the excess velocity. In a second step, we clarify the relevant influential parameters on these volume flows and their interconnection.

For the proposed modeling approach, we adapt a greybox model following Hangos and Cameron (2001). Our model is developed from engineering principles, hydrodynamic considerations (see Sec. 3) and well-defined equations, whereas the initialization of a part of the influential parameters is based on measured data. The underlying relations can be described as an intermediate concept between a black box (completely based on measurement data) and a white box model (based only on analytically well-known equations and engineering principles).

As depicted in the previous chapter (Sec. 3), we assume a droplet flowing through a rectangular microchannel with its properties: The thin wall film cross-section $A_{f}$ is determined by Eq. 12 and depends on the channel height $H$, the channel aspect ratio $ar$ and the film thickness $\delta$. To determine $\delta$, we apply the model of Han and Shikazono (2009), which holds for $Ca<0.2$.

The relative volume flow through the gutter ($Q_{g}^{rel}$) is derived from the pressure difference along the gutters as suggested by Abiev (2017). Therefore knowledge of the relevant pressures is crucial. It can be determined from the dynamic interface deformation caused by the moving liquids through the gutters in flow direction (Abate et al., 2012). Stagnant droplets have a static cap shape with a circular outline according to Musterd et al. (2015). When set into motion, the moving liquids exert forces onto the interface and cause a dynamic shape deformation (Luis A. M. Rocha and João M. Miranda and Joao B. L. M. Campos, 2017).

A model proposed by Mießner et al. (2019) allows to approximate the droplet shape and gutter diameter. The model implies that the deformation difference of the dynamic droplet cap shape between the droplet front and back results in a change of the gutter radius from the static shape. The cross-sectional gutter area $A_{g}$ widens asymmetrically from back to front with the growing gutter radii to accommodate the relative volume flow $Q^{rel}_{g}$ of the continuous phase through the gutter. The gutter entrance at the droplet front is therefore larger than the gutter exit of the droplet back. Utilizing their model, the dimensionless radius of these gutters can be calculated from the flow related curvature of the droplet.

The gutter radius $k_{g,i}$ is therefore defined as a fraction of the droplet height $h$. For the case of present wall films the term is expressed as follows:

[TABLE]

In order to simplify the geometry, we define a mean gutter radius $\overline{R_{g}}$ using Eq. 20:

[TABLE]

which is also used for the mean cross-sectional area $A_{g}$ derived with Eq. 11.

In previous work (Helmers et al., 2019), we introduce a quantification of the droplet cap deformation with an elliptic approximation of the cap outline

[TABLE]

[TABLE]

The ratio of the semi-major $a$ and semi-minor axis $b$ of the droplet cap curvature is introduced as deformation ratios $k_{c,i}$ at the droplet front and back. They become $k_{c,f}=k_{c,b}=1$ when describing the static circular droplet cap shape of the droplet front and back. For the dynamic cap shape, under the influence of the moving liquids, the droplet front appears elongated $k_{c,f}>1$ and the back cap is compressed in flow direction $k_{c,f}<1$. Hence we introduce a correlation to describe these relations: The cap curvature only depends on the $Ca$-number for moderate flows ($Re<5$):

[TABLE]

With the correlation approach from our recent publication and the model from Mießner et al. (2019), we are able to calculate the Laplace pressure at the gutters. We assume a linear connection between the gutter front and the back of the droplet body, since the curvature of the gutter in flow direction is negligibly small. In this case, the mean interface curvature at the gutter entrance and its exit depends on the gutter radii only and the flow induced Laplace pressure difference equals:

[TABLE]

[TABLE]

Those deformation related pressures at the droplet front and back provide a link to the driving pressure difference $\Delta p_{LP}$ along the gutter length in the flow direction. Due to symmetry, the static terms cancel out:

[TABLE]

Using the dimensionless expression from Eq. 20 this results in

[TABLE]

for the flow induced pressure difference as a driving force.

The relative volume flow rate through the four gutters $Q^{rel}_{g}$ can be modeled as a laminar pressure driven flow ((Bruus, 2008)) and it is linked to this pressure difference with a hydrodynamic resistance $\Omega$:

[TABLE]

The hydrodynamic resistance $\Omega$ is defined by Ransohoff and Radke (1988) and Shams et al. (2018) as

[TABLE]

Besides the mean gutter length $\overline{l_{g}}$, herein $\beta$ is a dimensionless factor, that represents the geometrical obstructions of the gutter flow, as well as the viscous coupling of both flow phases. In accordance with the simulation results of Shams et al. (2018), we declare an influence from the viscosity ratio $\lambda$ of the flow phases to take care of the viscous coupling effects:

[TABLE]

The mean gutter length $\overline{l_{g}}$ can be derived from the droplet length $l_{d}$, if the gutter distance $\Delta x_{g,i}$ from the caps is subtracted:

[TABLE]

The gutter distance from the front and back droplet tip $\Delta x_{g,i}$ was defined by Mießner et al. (2019) as

[TABLE]

The above stated considerations lead to our final description for the relative volume flow through the gutters

[TABLE]

Herein $\overline{l_{g}}$ and $\Delta p_{LP}$ depend on $\overline{R_{g}}$ as the preceding considerations show. Thus $\overline{R_{g}}$ has the most prominent influence feature besides $\beta$. Inserting Eq. 36 and Eq. 33 in Eq. 19 delivers

[TABLE]

expanding with $\frac{W}{W}$ and inserting $l_{g}$ (Eq. 34) we receive our final expression for the excess velocity:

[TABLE]

We point out, that our model is usable within the capillary regime $Ca<0.02$ since the model for the wall film-thickness from Han and Shikazono (2009), the analytic interface model from Mießner et al. (2019) and the droplet curvature correlation from our recent work (Helmers et al., 2019) are valid in this range. At higher $Ca$ in the viscous regime, different flow conditions with thicker wall films (Jose and Cubaud, 2014) as well as a higher influence of $Re$ are reported (Kreutzer et al., 2005).

5 Model Calibration

The final expression for the excess velocity (Eq. 38) depends on accessible data like $W$, $ar$ and $l_{d}$ as well as on $\overline{R_{g}}$, $\beta$, $k_{g,f}$, $k_{g,b}$, $\Delta x_{g,f}$ and $\Delta x_{g,b}$, whereas the latter are also related to the gutter radii $\overline{R_{g}}$. As shown in the previous section, the parameters can be calculated from the droplet cap curvatures $k_{c,f}$ and $k_{c,b}$ and via a measurement based model calibration. The 6 parameters $m_{cap,f}$, $m_{cap,b}$, $n_{cap,f}$, $n_{cap,b}$, $c_{cap,f}$, $c_{cap,b}$ influencing the cap deformation at the droplet front and back and additionally the dimensionless resistance factor $\beta$ with $m_{\beta},c_{\beta}$ and $n_{\beta}$ need to be adjusted.

For model calibration, we use the data-set presented in (Helmers et al., 2019) in combination with supplementary measurements in the data range of low $Ca$ and redefine $n_{cap}$ within an interval around 1. Nevertheless, for $Ca\rightarrow 0$ it represents the point of minimal surface energy and equates a spherical shape. This approach allows to improve the convergence of solvers. Out of further apriori considerations ($\beta>0$, $n_{c,f}\approx 1$), we additionally define the boundaries for the search space of the solver in Tab. 1. We allow the solver to adapt the correlation coefficient from our last work to the presented model, since correlations of measurement data unavoidably include measurement errors, that might bias the solver results.

Beneath the fixed boundaries of the search space, a hydrodynamic boundary condition is applied to improve the convergence of the used optimization algorithms. For a rising $Ca$, the difference between $k_{g,f}$ and $k_{g,b}$ must increase, because of the pressure difference between the droplet front and back increases with higher $Ca$ (Abiev, 2017) and the droplet front elongates, which leads to a larger front gutter radius, while the droplet rear flattens out. This is expressed by the gutter radius increase for an increasing $Ca$:

[TABLE]

The large number of influence parameters leads to a highly nonlinear optimization problem with numerous local minima. Thus most gradient-based algorithms are not suitable for this type of optimization problem, since they tend to converge to local optima. This would result in an enormous number of randomly initialized solver calls to cover the whole search space. Thus an optimization algorithm that is capable of global optimization e.g. stochastic and metaheuristic approaches are more favorable to cover the search space and solve the statistical part of the greybox model.

The quality of a solver result (e.g. deviation between measured data and estimation) is quantified by the loss-function of the problem. For the model calibration, we define

[TABLE]

with the weights of the individual properties $\omega_{1}$, $\omega_{2}$, $\omega_{3}$ following Tab. 2. Values with an upper index $\mathcal{M}$ denote values estimated by the model and values without an upper index represent the calibration data. The first two sums serve as calibration data-sets for the hydrodynamic flow properties, since they contain the flow related deformation and therefore the hydrodynamic influences. The differences of the velocity $u_{d}$ serves as a parameter for the actual droplet velocity and therefore the flow resistance $\beta$. This is necessary, since otherwise no representation for $\beta$ is available.

An additional factor $\psi$ is introduced to maintain the overall integrity of the model: The excess velocity depends on extrinsic measurable values such as the dimensionless quantities and geometrical properties as well as varying flow properties like the gutter length. Thus, it is appropriate to separate the measurable quantities from the model based quantities. In doing so, one can directly compare the quantities received from our correlation with measurement data adjusted for material and flow properties. The separation of those terms leads to the equilibrium function $\psi$ (Eq. 41 + 42). Herein Eq. 41 represents the data from our measurements and Eq. 42 only data from our modeling assumptions and geometry. For a well-adjusted model, the measured data for $\psi$ (upper equation) should correspond with the modeled values for $\psi$ (lower equation).

[TABLE]

The calibration procedure is based on a Genetic Algorithm (GA). Within the GA, every possible solution is emulated as a genetic code of population individuals. During the optimization process, the different solutions (individuals) can be combined (mated) to generate mixed solutions (children) with a combined genome. The decision, which solutions are actually combined, is based on the error (fitness value of loss function) of the solution. Within a so-called ranking sampling, the best solutions combine stronger than weak solutions (survival of the fittest), which leads to an improvement of the over-all population over the generations. Like in natural populations, random mutations of the genome can improve the overall fitness of a population. Transferred to a solver this means, that the population is able to leave local optima, if mutated individuals (solutions of partly random parameters) have a higher fitness value and therefore significantly change the genome pool of the population.

For good optimization results it is necessary to emulate a sufficiently sized genome pool, thus a high number of emulated individuals is preferred. This in turn results in a massively increased calculation demand, because for each iteration step every single individual and the children must be evaluated (Whitley, 1994). Additionally, the genetic algorithm is typically performed several times to identify local optima.

To decrease the amount of time-consuming iteration steps of the GA and thereby reduce overall calculation time, we utilize a three-step stochastic and gradient free approach:

In the first step, a random population at the feasible borders of the problem is generated for the Genetic Algorithm and a genetic optimization performed. The convergence point of the GA is initialized via Latin Hypercube Sampling processed by a following fast-converging Pattern Search Algorithm (PSA) (Davey, 2008), that results in an improved minimum as the final convergence point. The properties of both algorithms are shown in Tab. 3. The algorithm finally merges at the values shown at the end of this section (Tab. 4). The results are discussed in the following.

For the flow induced cap curvature, we again find the corresponding interrelation from our last work (Helmers et al., 2019). For both cap deformation ratios $k_{f}$ and $k_{b}$ an exponential behavior as a function of the $Ca$-number is visible (Fig. 4). This proves the assumption of our previous work (Helmers et al., 2017) and agrees to Mießner et al. (2019). Our additional boundary condition (Eq. 39) is satisfied for all values as the graphs for the gutter radii $k_{g,i}$ show.

Rising viscous forces, indicated by a rising $Ca$-number, deform the droplet interface more strongly. At static conditions ($Ca\rightarrow 0$) the cap curvatures are in correspondence with Musterd et al. (2015) roughly circular ($k_{c,f}=k_{c,b}\approx 1$). The back cap is compressed with respect to the main flow direction ($k_{c,b}<1$) and the front cap elongated ($k_{c,f}>1$) if viscous forces rise in comparison to the static case.

The size of the gutter radii rises with the increasing influence of the viscous forces, since the bypass flow in the gutter increases and needs to be accommodated by the gutters. The gutter entrance is always larger in diameter than the exit radius ($k_{g,f}>k_{g,b}$).

For the resistance factor $\beta$ no fitting data is available since it can not be measured directly. Therefore $\beta$ is fitted based on the velocity data from high-speed camera measurements and the application of parameters for the droplet deformation. The resulting droplet velocities in comparison with the measured values are shown in Fig. 5. The measurements fit reasonably well within the range of inevitable velocity fluctuations caused by Taylor flow stability of +/- 10 % range described by (Vansteijn et al., 2008).

All adaption coefficients for the proposed model are summarized in Tab. 4.

6 Model Validation

A a first validation step the calibration functions $\psi$ and $\psi^{\mathcal{M}}$ are considered. For a hydrodynamic well-adjusted model, both function should coincide and our model based shape deviations ($\psi^{\mathcal{M}}$) equals the combination of measured properties ($\psi$). The corresponding data and the values for our model agree well at high $Ca$ numbers (Fig. 6), whereas a deviation for lower $Ca$-numbers can be observed. This can be explained by the fact, that the excess velocity itself is a relative quantity and it is therefore stronger influenced at lower absolute values (low $Ca$-number). Thus an inevitable constant measurement deviation for velocity and volume flows caused by the experimental equipment results in a higher error for low excess velocities. Especially for $Ca<10^{-4}$ the resulting volume flows are situated at $Q_{0}\approx 2\cdot 10^{-6}l/min$ and even minor deviations lead to high errors for the excess velocity. Thus, we consider our model approach to represent the measurements reasonably well and our fit coefficients to be valid.

Besides the hydrodynamic validation, our assumptions for $\beta$ are compared with available simulation data. We find a dependence of the gutter flow resistance $\beta$ on the viscosity ratio $\lambda$. The latter can be interpreted as an indicator for the viscous coupling of both flow phases (Fig. 7). In case of a highly viscous continuous phase ($\lambda<1$) the strongest velocity gradients are found inside the droplet, while for a viscous disperse phase ($\lambda>1$) larger velocity gradients and therefore energy dissipation is found inside the gutter-flow, resulting in a larger $\beta$.

This approach agrees with the simulations from Shams et al. (2018), who improved the model of Ransohoff and Radke (1988) by introducing the viscous coupling of the disperse and the continuous phase. For our case of $\lambda=0.1-1.4$ and a contact angle $\theta=0$ Shams et al. (2018) report a $\beta$ between 20-30 for a co-current flow. Our values are shifted by a constant offset, while the slope and therefore the dependence on $\lambda$ is very similar. We consider this to be caused by the use of a different flow field specifications. Shams et al. (2018) describe a concurrent flow in Eulerian specification, while in this work we determine $Q_{g}^{rel}$ within a Lagrangian flow specification. The coordinate transformation thus can only change the offset of the function, while the hydrodynamic influence (the slope) must remain identical. Additionally within the simulation of (Shams et al., 2018), they assume a contact line between disperse phase, continuous phase and the wall in their problem definition. Although for the solution shown in Fig. 7 the contact angle for the continuous phase is nearly 180 °, the existence of a contact line introduces an additional resistance. Thus regarding $\beta$ we consider our model valid.

The integrity of the model itself and the correlation for $\beta$ has been successfully proven, but due to the mentioned experimental restrictions, we can not directly compare the modeled and experimental determined excess velocities $u_{ex}$. Instead we compare the results of our model to different published approaches.

A suitable correlation for the prediction of the excess velocity in the first regime was introduced by Jose and Cubaud (2014), who identify the $Ca$-number and the ratio of the droplet length $\frac{l_{d}}{H}$ as characteristic properties (Fig. 8). It has to be mentioned, that the model of Jose and Cubaud (2014) for non-wetting droplets ends at $\frac{l_{d}}{HCa^{-1}}>600$ , since for larger values they observe the disperse phase to wet the channel walls. This results in intensified dissipation and a higher pressure drop and thereby inhibits a gutter flow from the droplet front to the back. This phenomenologically equals the second flow regime as mentioned by Wong et al. (2004), but lacks the thin wall film and results in a much larger pressure drop. Furthermore, the viscosity ratio $\lambda$ is not included in their correlation.

A comparison of our deterministic model with Jose and Cubaud (2014) correlation Fig. 8 shows good agreement. At very low $Ca$-numbers ($Ca<10^{-4}$), our model results in a slightly increased excess velocity. We regard this behavior of the model as not physical. The effect results from the mathematical counterplay of the terms $\lim\limits_{Ca\rightarrow 0}{\frac{1}{Ca}}=\infty\leftrightarrow\lim\limits_{Ca\rightarrow 0}k_{g,f}-k_{g,b}=0$ within Eq. 38. In order to achieve a stable solver convergence, we accepted a small residual deviation for the static case $Ca\rightarrow 0$ for the front and back shape of $0.1\%$. Due to the relative character of the excess velocity, this unfolds a significant influence at low $Ca$ values. Unfortunately, additional measurement validation concerning the shape deviation at very low $Ca$ is not possible in our experimental design, since the expected shape deviation is smaller than the blurriness of the interfacial area in the images itself and therefore lies within the measurement deviation.

7 Discussion

The interfacial area of a Taylor droplet in rectangular channels can be divided into the front and back cap regions, the wall films and the gutter interface. Neglecting the caps, the main momentum input into a droplet is transferred across the wall film and the gutter interface area. An increasing channel aspect ratio and droplet length result a growing wall film area, i.e. an enlarged dissipation interface.

As we showed in Sec. 2, the behavior of the Taylor droplet’s excess velocity can be parted in two possible regimes and the viscosity ratio $\lambda$ has a strong influence on the hydrodynamic mechanisms (Sec. 2).

Within the first regime the fluid in the gutter flows slower than the droplet, exerting a drag force and leading to a positive excess velocity $u_{ex}>0$. These drag forces influence the droplet shape, leading to a flattened droplet back and elongated droplet front. We characterize this shape variation with a correlation (Sec. 5). As it can be seen from the measurements, our data falls into this flow regime, as our resulting droplet shape indicates in accordance with Wong et al. (1995).

The influence of the droplet length in the plug flow regime and $\lambda<1$, where larger droplets have an $u_{ex}\approx 0$ is in agreement with Jakiela et al. (2011) and their later publication Jakiela et al. (2012). Additionally, they find an elongation of the dynamic droplet length in comparison to the static droplet length, which is also covered by our shape correlation, since for rising the $Ca$-number the droplet front elongates stronger than the droplet back is compressed. Recent published simulations by Kumari et al. (2019) show, that also for larger $Re$-numbers $u_{ex}\approx 0$.

Our concept of deriving the excess velocity from the gutter pressure drop, that is inverted to the flow direction (larger pressure at the front gutter entrance) is in accordance with Abiev (2017). Nevertheless, the averaged pressure is still higher at the droplet back than on the droplet front due to the overall droplet pressure drop, since the droplet needs a driving force for its translation.

For the second flow regime identified by Wong et al. (1995), where the viscous dissipation in the film and droplet leads to a bypass flow from the droplet back to the front and therefore a negative droplet excess velocity (for large $\lambda$ and long droplets), our model can be adapted, if the gutter-shape-difference term is revised or a resistance coefficient for the film is added. As the viscosity ratio $\lambda$ rises, more momentum will be dissipated via the wall films. This extra momentum is dissipated at the gutter interface, which results in a slower droplet velocity, forcing the continuous phase to bypass the droplet reversely. Since we are unable to establish a water-in-oil two-phase flow for a $\lambda\ll 1$ in our experimental setup due to the hydrophilized channel walls, droplet shape correlations for the case of $\lambda\gg 1$ should be performed in future work. Therefore, although we assume a systematic inversion of the gutter radii ratio back to front as a consequence of the reversed gutter flow direction, we like to mention that with the current shape correlation our model only works for the case of low viscous disperse phase ($\lambda\lessapprox 1$) like gas/liquid or low viscous oil/water flows.

The influence of the channel aspect ratio as mentioned by Wong et al. (2004) is incorporated in our model: a higher aspect ratio results in lower excess velocities. This can be explained by the larger drag forces acting on a larger relative wall film area caused by the flattened channel geometry.

The comparison with the most prominent approaches shows, that our model and the chosen influential parameters are valid for moderate and small viscosity ratios. The excess velocity is determined by viscous dissipation within the droplet and the gutters, as well as the drag of the thin wall films. The relation is characterized by the $Ca$-number, viscosity-ratio $\lambda$, the dimensionless gutter-length $l_{g}$, the aspect ratio $ar$ and the wall-film thickness $A_{f}$. Furthermore, the proposed model can close the gap for $\frac{l_{d}}{HCa^{-1}}>600$ and allows the calculation of the excess velocity for moderate $Ca$-numbers ($Ca<0.02$).

8 Conclusion

In this work, we developed a model to determine the relative droplet velocity of Taylor flows in square microchannels for moderate $Ca$-numbers and low to moderate viscosity disperse phase ($\lambda\lessapprox 1$). We base our model on the relative volume flow through the gutters, as well as the wall film thickness. The flow through the gutters is determined from the pressure drop described by the Laplace-Pressure difference between the gutter entrances.

Our model uses the gutter radii to obtain the resulting pressure gradient that drives the continuous phase through the gutters. We use measurements at different $Ca$ and $Re$ in surfactant-free fluid system from a previous publication (Helmers et al., 2019) to derive the radii at the gutter entrances from the surface shape model proposed by Mießner et al. (2019) and calibrate the model parameters.

Our model is successfully validated with an intrinsic approach comparing the congruence of measurement data and calibrated model parameters. Additionally, we successfully compared our model to the phenomenological correlation of Jose and Cubaud (2014).

For the future, our model should also be validated for aspect ratio differing from unity, since the influence of aspect ratio has been integrated into our model, but has not been validated so far. Additionally, the influence of surfactants and highly viscous droplets ($\lambda\gg 1$) on the excess velocity should be investigated to extend the model, since an excess velocity $u_{ex}<1$ was not included into the model so far. We suggest to include this function in the modelling of the wall film resistance. Especially local spatially resolved measurement techniques, e.g. µ-PIV measurements should be appropriate for this task.

\authorcontributions

conceptualization, T.H. and U.M.; methodology, T.H. and U.M.; software, T.H., P.K. and U.M.; validation, T.H.; formal analysis, T.H. and U.M.; investigation, T.H.; resources, T.H.; data curation, T.H.; writing—original draft preparation, T.H.; writing—review and editing, U.M. and J.T.; visualization, T.H. and P.K.; supervision, U.M. and J.T.

\funding

This research received no external funding.

\conflictsofinterest

The authors declare no conflict of interest.

\abbreviations

The following abbreviations are used in this manuscript:

\appendixtitles

yes \appendixsectionsone

Appendix A Considerations for $u_{ex}$

Rearranging Eq. 18 leads to the equation

[TABLE]

Our measurements show in agreement with Jose and Cubaud (2014) excess velocities with values $u_{ex}<0.4$ for $Ca<0.2$. Additionally we can assume $\frac{A_{f}}{A_{ch}}<0.005$ as shown in Fig. 9. Thus one can say $u_{ex}\frac{A_{f}}{A_{ch}}<0.002$ and therefore it can be considered small of higher order and be neglected.

\externalbibliography

yes