# Slug Self-Propulsion in a Capillary Tube Mathematical Modeling and Numerical Simulation

Authors: M. I. Khodabocus; M. Sellier; V. Nock
Journal: Advances in Mathematical Physics
 (2016)

---

## Abstract

A composite droplet made of two miscible fluids in a narrow tube generally moves under the action of capillarity until complete mixture is attained. This physical situation is analysed here on a combined theoretical and numerical analysis. The mathematical framework consists of the two-phase flow phase-field equation set, an advection-diffusion chemical concentration equation, and closure relationships relating the surface tensions to the chemical concentration. The numerical framework is composed of the COMSOL Laminar two-phase flow phase-field method coupled with an advection-diffusion chemical concentration equation. Through transient studies, we show that the penetrating length of the bidroplet system into the capillary tube is linear at early-time regime and exponential at late-time regime. Through parametric studies, we show that the rate of penetration of the bidroplet system into the capillary tube is proportional to a time-dependent exponential function. We also show that this speed obeys the Poiseuille law at the early-time regime. A series of position, speed-versus-property graphs are included to support the analysis. Finally, the overall results are contrasted with available experimental data, grouped together to settle a general mathematical description of the phenomenon, and explained and concluded on this basis.


---

## Body

## 1. Introduction

When two physically distinct but miscible liquid droplets are placed in a capillary tube, a key feature that is commonly observed is the bulk motion of the composite droplet under the action of capillarity until complete mixing is attained [1]. In Figure 1 is given a freehand sketch of this scenario. This is an example of motion due to (chemical)Marangoni effect [2–4]. The motion of the bidroplet system is caused by local variations of its advancing and receding surface tensions that are in turn caused by differences in its chemical compositions. Out of this simple scenario grew two kinematical concepts, namely, thepenetrating length and thepenetrating rate. When both concepts are grouped together, we attain, by inference,slug self-propulsion mechanism—a juxtaposition of two miscible liquid droplets. An important facet of this driving mechanism is that it is simple and so could be delicately adjusted to manipulate liquid droplets in small-scale fluidic devices. Earlier scientists have studied the phenomenon experimentally and have concluded that it might be used as a means of self-propelling, transporting microscale substances at a respectable penetrating length and rate, respectively, of orders of centimeters, centimeters per second [1, 5, 6].Figure 1
Illustration of the slug self-propulsion phenomenon in a capillary tube.Bouasse [7] inquired experimentally on capillarity and interfacial phenomena. He found that a droplet predeposited in a conic tube undergoes motion towards the narrow region because of the Laplacian pressure difference induced on both its menisci. Weislogel [8] investigated the dynamics of a liquid droplet in a chemically asymmetric tube, by treating the tube to be successively hydrophilic and hydrophobic. He found the droplet to move towards the hydrophilic side. Dos Santos and Ondarçuhu [9] argued that a droplet containing a reactive substance remained in its state of motion as long as such reactive agent is operative—the dynamics of reactive wetting on solid surfaces were discussed by De Gennes [10].Historically, it may be that the first documented study goes back to 1712, when Hauksbee [11] investigated the behaviour ofa drop of oil of oranges sandwiched in between two nonparallel glass plates. In that paper, two points obviously follow in the light of this phenomenon. First, the fact that the motion of the drop reversed when the slope of those two plates was reversed led Hauksbee to believe that the direction of motion is towards the narrow region. Secondly, the fact that his experiment was reproducible in vacuum further informed Hauksbee that such motion is solely due to interfacial phenomenon.Currently, existing (active) techniques of manipulating discrete droplets includeelectrowetting [12–17],dielectrophoresis [18],thermocapillary [19–22],surface acoustic wave (SAW) [23, 24], andpressure-driven based-channel systems [25]. Thanks to these sophisticated techniques, the fields of nano- and microfluidics have seen over the past decade considerable development.Self-propulsion phenomena as well as coalescence phenomena are found to have generated considerable interest and many features have been correctly described and explained as well. Some fundamental questions were addressed by De Gennes [10, 26, 27] and Shanahan and De Gennes [28]. Shanahan and De Gennes stressed that the motion of the system is switched on mainly due to the slight perturbation communicated to the system at the time when the second droplet is brought into close contact to the predeposited one. Furthermore, the authors pointed out that for motion to occur the system must necessarily be asymmetric and that such motion is maintained only insofar that such character is observed in the system.Recently, there has been a growing interest in exploiting droplets self-propulsion phenomenon, thereby avoiding the arduous fabrication procedures required [6, 29, 30]. For instance, Karpitschka and Riegler [31] gave an account on the coalescence of liquid droplets in two respective regimes, namely, thefast coalescence ordelayed coalescence. In the paper of Karpitschka and Riegler, the authors pointed out that the coalescence phenomenon is predominantly governed by surface tension, regardless of whether it is fast or delayed coalescence. For instance, when the difference δγ in the surface tensions of the working droplets is such that δγ≈3 mN·m−1, delayed coalescence phenomenon is observed from the bidroplet system and, for all other cases such that δγ<3 mN·m−1, fast coalescence phenomenon is observed. In a subsequent paper, Karpitschka and Riegler [32] explore the study of delayed droplets coalescence phenomenon on a combined experimental and theoretical basis.On the one hand, the self-propulsion of liquid droplets has also been addressed by Sinz and Darhuber [6] and Sellier et al. [29, 30]. In the works of Sellier et al., the authors examined the self-propulsion mechanisms in unconfined systems, that is, sessile droplets configuration. The principles of manipulating discrete droplets using active methods were considered by Darhuber, Valentino, Davis, and Troian [33, 34]. In the paper of Sinz and Darhuber, the authors addressed the self-propelling of surfactant droplets in chemically confined systems. Therein, Sinz and Darhuber argued that this passive driving mechanism might be used to transport microscale substances.Lunati and Or [35] undertook the problem of gravity-driven slug dynamics in capillary tubes. The authors have answered questions regarding the origin, strength, and effects of the drag-forces developed at the contact lines. They also studied the kinematics of the slug, contrasting the numerical simulation results with the experiment results of Bico and Quéré [36] on falling of viscous slugs in vertical capillary tubes. Bico and Quéré deduced a set of scaling relationships characterising the steady motion of the slug. On the experimental side, the earliest laboratory experiments wherein the subject matter based on purely confined systems has been explored thoroughly were those of Bico and Quéré [1, 5]. These are the most important works in close connection with the subject at hand. In these, more work was endeavoured on the dynamics of trains of droplets in horizontal capillary tubes in both wetting and drying conditions. For instance, in the former paper of Bico and Quéré [5], which constitute their preliminary ground-based experiments, they argued mainly on the coating properties of trains of juxtaposed liquid droplets in terms of geometrical and viscous effects—another very interesting field of inquiries; they moreover drew attention on their versatility towards integrated microfluidic devices. For example, Bico and Quéré [1] have studied the response of the system in more complex geometries and displayed more features on a combined experimental and theoretical basis.Since then, only few scientists have begun to study the subject on theoretical and numerical sides. Furthermore, studies numerically based on the slug self-propulsion phenomenon do not appear to have received much attention. To sum up, it may be concluded that we found no scientist has come to ponder on the self-propelling phenomenon in confined system quantitatively and to arrive at a more rational conception of it.We have studied carefully the dynamics of a two miscible fluids in a capillary tube and, accordingly, found the law describing the position (orpenetrating length) of the slug into the capillary tube is well described by the following relationship: Xst≈η11-exp-η2t-t0, where η1, η2>0 are property-dependent constants. Furthermore, it is observed at the early-time that the slug velocity (orrate of penetration),that is, Us=∂tXst, obeys aPoiseuille law: Us∝Rγ1,0-γ2,0/4μls. In this expression, R designates the radius of the capillary tube; γ1,0, γ2,0>0 are the surface tensions of the miscible fluids; and μ,ls>0 are the dynamic viscosity and length of the slug. We came to these laws numerically by collecting data on the slug position as a function time and the maximum travelled distance. An inspection of the form of the property Us (discussed hereinafter) evidences at once that small and low viscous substances are speedier than large and highly viscous ones.In this paper, an account of the mathematical development is given. The paper is organised as follows: Section2 is devoted to mathematical formulations and Section 3 to the numerical implementation of the mathematical framework. Results from numerical experiments and analysis are presented in Section 4: the motions of the slug with respect to time and property constants are discussed in Sections 6 and 7, respectively; the results are contrasted with the papers of Bico and Quéré [1, 5]. In Section 8, a power-law model describing the position of the slug with respect to time and the leading property constants involved is established. The overall results are finally argued in Section 9.


## 2. Mathematical Formulations

### 2.1. Preliminaries

First, some aspects, definitions, and notations are in order. Truly speaking, we examine here athree-fluid flow problem, the properties of which are such that two of which are miscible phases but immiscible with the third one. In this context, their respective domains of definition are written as (1)Ωv-∪Ωv+=defr,z:0≤r≤R,0≤z<zΓ2∪zΓ1<z≤L,Ω1∪Ω2=defr,z:0≤r≤R,zΓ2≤z≤zΓ1,where r=r,zT is a point of either the set Ωv-∪Ωv+ or Ω1∪Ω2. The scalar fields zΓit (with i=1,2) are the respective time-dependent positions of the (front and rear) interfaces of the slug. Obviously, the sets Γ1=Ω1∪Ωv+ and Γ2=Ω1∪Ωv- denote those time-dependent interfaces. These physical domains, sets (1), are clearly illustrated in Figure 2. The choice of Figure 2 is motivated by that of Figure 1. In view of the foregoing definitions, it thus appears evident to regard Ω1∪Ω2⊂R2 merely as atwo-miscible-liquid region and Ωv±=Ωv-∪Ωv+⊂R2 the (disjoint) union of its rear and frontvapour phase regions. In fact, with these underlying definitions in mind, both phases can be related without loss of generality as thus (2)volΩv-∪Ω1∪Ω2∪Ωv+=volΩ1∪Ω2+volΩv-∪Ωv+,where vol· denotes thevolume of its arguments. Furthermore, we characterise by the scalar fields ϕ12,ϕv∈0,1 the functionals describing thevolume fraction of fluids in Ω1∪Ω2⊂R2 and Ωv-∪Ωv+⊂R2, respectively. Consequently, by virtue of the above postulation, (2), it results in(3)1=volΩ1∪Ω2volΩv-∪Ω1∪Ω2∪Ωv++volΩv-∪Ωv+volΩv-∪Ω1∪Ω2∪Ωv+=ϕ12+ϕv,where ϕ12Γi≃ϕvΓi≃1/2 (with i=1,2). In regard to these facts, the problem under consideration now reads anImmiscible Two-Phase Flow problem. At this point, we have made the essential geometrical features prominent and, in what follows, a two-phase flow concentration-dependent surface tension problem of a viscous liquid slug in a capillary tube is undertaken.Figure 2
Illustration of a slug composed of two miscible fluidsΩ1∪Ω2, both of which bounded partly by the capillary tube Γs and by Γk (with k=1,2) by a still vapour phase Ωv±. The lengths L, l1, l2>0 denote the length of the capillary tube and droplets 1 and 2. Bold arrow denotes the direction of its motion.


### 2.2. Transport of Concentration by Diffusion and Convection

Denote by the scalar fieldClr,t the chemical concentration of a (highly) volatile substance filling the region Ω2 initially (note that the other substance filling the region Ω1 is taken to be lowly volatile). To establish the transport equation for the chemical concentration function Clr,t in Ω1∪Ω2, we reason in the following way.Initially, postulate the concentration within the slug to be such that(4)Clr,t=i-1inΩi×t=0i=1,2;and presuppose anormally imposed flux φi≥0 on the corresponding boundary Γi  (with i=1,2,s) by its external surroundings. Then, by virtue of these assertions, Clr,t is simply described by the equation of chemical species transported by diffusion and convection (CD-model, say), as(5)∂tCl+divClu-divD21∇Cl=0in⋃i=1,2Ωi×0,T;-ni·Clu-D21∇Cl=φionΓi×0,Ti=1,2,s;Clr,t=i-1inΩi×0i=1,2;for given T>0. In (5), the proportionality factor D21>0 is thediffusivity of fluid 1 into 2; the vector field nir,t denotes the outward normal unit vector corresponding to the boundary Γi (with i=1,2,s). It may be pointed out in passing that this convection-diffusion model,that is, (5), merits further attention, for it should also agree with the subsequentchanges in position of the slug. This point is clarified in a later section.


### 2.3. Inclusion of Concentration-Dependent Surface Tension

We now discuss proper forms of interfacial forces at the boundariesΓ1, Γ2. Supposing a solution for the scalar field Clr,t is granted, we can then define γi(CΓil) as thegrowth, i=1 (resp.,decay, i=2) of interfacial effects in terms of CΓil=ClΓi at their respective boundaries Γ1, Γ2. In order to deduce such proper forms for γi(CΓil), we shall in the first place consider the following facts: initially, the domains Ω1, Ω2 are characterised to be such that Cl=0 in Ω1 and Cl=1 in Ω2, and, finally, both are expected to be in equilibrium (chemically), attaining the value C∞l∈0,1, say; thence, the correct conditions to be imposed upon CΓklr,t (with i=1, 2) simply write as(6)CΓil=i-1onΓi×t=0i=1,2;CΓil=C∞lon⋃i=1,2Γi×t⟶∞.Taking (6) into consideration along with an additional assertion that the surface tensions evolve linearly, one arrives at(7)γiCΓil=C∞l-CΓilC∞l-i-1γi,0+CΓil-i-1C∞l-i-1γi,∞onΓi×0,Ti=1,2.The linearity of γi(CΓil) (with i=1,2) on the scalar field Clr,t is here motivated by the works of Carles and Cazabat [37, 38] on the spreading dynamics of Newtonian droplets. Therein, Carles and Cazabat [37] found that the surface tension of those mixed substances varied roughly linearly with the chemical concentration of the mixture. Equation (7) is physically reasonable, for Cl∈[0,C∞l] in Ω1 and Cl∈[C∞l,1] in Ω2. However, out of (7) grew one important question which is somehow unclear.Question 1.
What are the magnitudes ofC∞l, γi,∞>0 (with i=1,2)?The answer of this question is found in the following line of reasoning. First, it is evident that(8)liml2→0⁡C∞l=0,liml2→∞⁡C∞l=1for  everyl1=Const.,t∈0,T,where the geometrical property l1>0 stands for the size of Ω1 along the z-axis (see Figure 2). That is, the (chemical) equilibrium state is somewhat dependent upon the geometrical properties of the slug, as found in the paper of Bico and Quéré [1]. Thus, the expression which would match these conditions is C∞l=l2/l1+l2, since Cl=i-1 in Ωi (with i=1,2) at time t=0. Note that C∞l=1/2 is also an admissible value but only for a particular case so to speak when l1=l2, for it neglects the geometrical property of the slug, particularly, when l2→0, ∞.In regard to these arguments, one then concludes that the lesser the lengthl2>0 the higher the concentration C∞l∈0,1. It now remains to settle proper forms for γi,∞>0 (with i=1,2). First, we note that(9)limt→∞⁡l1γ1+l2γ2l1+l2=l1γ1,∞+l2γ2,∞l1+l2=l1γ1,0+l2γ2,0l1+l2on⋃i=1,2Γi;γ2-i,∞≡γ2-i,0whenl1+i=0i=0,1.Hence, these suggest to consider abarycentric-based model, namely, of the form (10)γi,∞=l1γ1,0+l2γ2,0l1+l2,onΓii=1,2,because γ1,∞=γ2,∞ (equilibrium surface tension) when t→∞. Clearly, there is no loss of generality in doing so, for not only does (10) agree well but it is also a physically acceptable expression. We are now well prepared to set forth the governing equations; this is the purpose of the next section.


### 2.4. Two-Phase Flow Phase-Field Model

In Sections2.2 and 2.3, we covered much ground on diffusive and interfacial effects. The present section takes those effects into account and models the slug self-propulsion mechanism using thephase-field approach. In a sense, we have beforehand the motion of a two-miscible-liquid system partly in contact with an immiscible one, the vapour phase region. Therefore, the starting point is naturally based on the classical Navier-Stokes equations (NS-model, say) for incompressible flow with interfacial forces coupled with the so-calledCahn-Hilliard evolution equation (CH-model, say). In the absence of gravitational field, those equations in vector invariant forms read as follows: (11)∂tρu+divρu⊗u-divTp,u-G∇ϕ=0in⋃i=1,2,vΩi×0,T,∂tϕ+divϕu-divϖ∇G=0in⋃i=1,2,vΩi×0,T,divu=0in⋃i=1,2,vΩi×0,T,where the stress tensor Tp,u, constructed with the pressure and velocity fields pr,t, ur,t reads as follows:(12)Tp,u=-pI+2μSu,whereSu=12∇u+∇uT.It is understood that Ωv=Ωv-∪Ωv+ (=Ωv±). In (11) and (12), the property constants ρ,μ>0 characterise themixed density and dynamic viscosity of the slug; both are defined in terms of the material properties of the two phases, and the smooth (phase-field) function ϕ12r,t∈0,1 is β=βv+β2-βvϕ12 (with β=ρ, μ). The function Gϕ defines the so-calledbulk chemical potential as(13)Gϕ=-ζε2divε2∇ϕ-ϕ2-1ϕ,whereϕr,t∈-1,1.In (13), ε>0[m] is acapillary width that scales with the thickness of the boundaries Γ1, Γ2; ϖ>0 [m3·s·kg−1], themobility term which determines the time scale of the Cahn-Hilliard diffusion; and ζ>0[N], themixing energy density. For more details see, for instance, the papers of Yue et al. [39] and Jacqmin [40] as well. We draw attention to the fact that ϕ12≃1+ϕ/2, and, therefore(14)ϕ12=min⁡max⁡1+ϕ2,0,1,whereϕr,t∈-1,1.In (14), min· and max· are operators which give the minimum and maximum of their respective arguments. Throughout the following discussion, the phase-field function ϕr,t∈-1,1 is used instead of ϕ12r,t∈0,1. It now remains to supplement the appropriate boundary conditions on Γ1,Γ2,Γs. At theinlet andoutlet boundaries azero pressure is imposed;that is, the system containing the domains Ω1, Ω2 is embedded in astill environment Ωv.Wetted wall boundary conditions are imposed along Γs; the contact angles are such that θ1≡θ2 (≡π/3 [rad]). The choice θk≡π/3 [rad] (with k=1,2) is justified from the fact that the advancing and receding boundaries of the miscible bidroplet system Ω1∪Ω2 are convex curves, following the works of Bico and Quéré [1, 5]. More features based partly on the field ϕr,t∈-1,1 are discussed in Section 3.


## 2.1. Preliminaries

First, some aspects, definitions, and notations are in order. Truly speaking, we examine here athree-fluid flow problem, the properties of which are such that two of which are miscible phases but immiscible with the third one. In this context, their respective domains of definition are written as (1)Ωv-∪Ωv+=defr,z:0≤r≤R,0≤z<zΓ2∪zΓ1<z≤L,Ω1∪Ω2=defr,z:0≤r≤R,zΓ2≤z≤zΓ1,where r=r,zT is a point of either the set Ωv-∪Ωv+ or Ω1∪Ω2. The scalar fields zΓit (with i=1,2) are the respective time-dependent positions of the (front and rear) interfaces of the slug. Obviously, the sets Γ1=Ω1∪Ωv+ and Γ2=Ω1∪Ωv- denote those time-dependent interfaces. These physical domains, sets (1), are clearly illustrated in Figure 2. The choice of Figure 2 is motivated by that of Figure 1. In view of the foregoing definitions, it thus appears evident to regard Ω1∪Ω2⊂R2 merely as atwo-miscible-liquid region and Ωv±=Ωv-∪Ωv+⊂R2 the (disjoint) union of its rear and frontvapour phase regions. In fact, with these underlying definitions in mind, both phases can be related without loss of generality as thus (2)volΩv-∪Ω1∪Ω2∪Ωv+=volΩ1∪Ω2+volΩv-∪Ωv+,where vol· denotes thevolume of its arguments. Furthermore, we characterise by the scalar fields ϕ12,ϕv∈0,1 the functionals describing thevolume fraction of fluids in Ω1∪Ω2⊂R2 and Ωv-∪Ωv+⊂R2, respectively. Consequently, by virtue of the above postulation, (2), it results in(3)1=volΩ1∪Ω2volΩv-∪Ω1∪Ω2∪Ωv++volΩv-∪Ωv+volΩv-∪Ω1∪Ω2∪Ωv+=ϕ12+ϕv,where ϕ12Γi≃ϕvΓi≃1/2 (with i=1,2). In regard to these facts, the problem under consideration now reads anImmiscible Two-Phase Flow problem. At this point, we have made the essential geometrical features prominent and, in what follows, a two-phase flow concentration-dependent surface tension problem of a viscous liquid slug in a capillary tube is undertaken.Figure 2
Illustration of a slug composed of two miscible fluidsΩ1∪Ω2, both of which bounded partly by the capillary tube Γs and by Γk (with k=1,2) by a still vapour phase Ωv±. The lengths L, l1, l2>0 denote the length of the capillary tube and droplets 1 and 2. Bold arrow denotes the direction of its motion.


## 2.2. Transport of Concentration by Diffusion and Convection

Denote by the scalar fieldClr,t the chemical concentration of a (highly) volatile substance filling the region Ω2 initially (note that the other substance filling the region Ω1 is taken to be lowly volatile). To establish the transport equation for the chemical concentration function Clr,t in Ω1∪Ω2, we reason in the following way.Initially, postulate the concentration within the slug to be such that(4)Clr,t=i-1inΩi×t=0i=1,2;and presuppose anormally imposed flux φi≥0 on the corresponding boundary Γi  (with i=1,2,s) by its external surroundings. Then, by virtue of these assertions, Clr,t is simply described by the equation of chemical species transported by diffusion and convection (CD-model, say), as(5)∂tCl+divClu-divD21∇Cl=0in⋃i=1,2Ωi×0,T;-ni·Clu-D21∇Cl=φionΓi×0,Ti=1,2,s;Clr,t=i-1inΩi×0i=1,2;for given T>0. In (5), the proportionality factor D21>0 is thediffusivity of fluid 1 into 2; the vector field nir,t denotes the outward normal unit vector corresponding to the boundary Γi (with i=1,2,s). It may be pointed out in passing that this convection-diffusion model,that is, (5), merits further attention, for it should also agree with the subsequentchanges in position of the slug. This point is clarified in a later section.


## 2.3. Inclusion of Concentration-Dependent Surface Tension

We now discuss proper forms of interfacial forces at the boundariesΓ1, Γ2. Supposing a solution for the scalar field Clr,t is granted, we can then define γi(CΓil) as thegrowth, i=1 (resp.,decay, i=2) of interfacial effects in terms of CΓil=ClΓi at their respective boundaries Γ1, Γ2. In order to deduce such proper forms for γi(CΓil), we shall in the first place consider the following facts: initially, the domains Ω1, Ω2 are characterised to be such that Cl=0 in Ω1 and Cl=1 in Ω2, and, finally, both are expected to be in equilibrium (chemically), attaining the value C∞l∈0,1, say; thence, the correct conditions to be imposed upon CΓklr,t (with i=1, 2) simply write as(6)CΓil=i-1onΓi×t=0i=1,2;CΓil=C∞lon⋃i=1,2Γi×t⟶∞.Taking (6) into consideration along with an additional assertion that the surface tensions evolve linearly, one arrives at(7)γiCΓil=C∞l-CΓilC∞l-i-1γi,0+CΓil-i-1C∞l-i-1γi,∞onΓi×0,Ti=1,2.The linearity of γi(CΓil) (with i=1,2) on the scalar field Clr,t is here motivated by the works of Carles and Cazabat [37, 38] on the spreading dynamics of Newtonian droplets. Therein, Carles and Cazabat [37] found that the surface tension of those mixed substances varied roughly linearly with the chemical concentration of the mixture. Equation (7) is physically reasonable, for Cl∈[0,C∞l] in Ω1 and Cl∈[C∞l,1] in Ω2. However, out of (7) grew one important question which is somehow unclear.Question 1.
What are the magnitudes ofC∞l, γi,∞>0 (with i=1,2)?The answer of this question is found in the following line of reasoning. First, it is evident that(8)liml2→0⁡C∞l=0,liml2→∞⁡C∞l=1for  everyl1=Const.,t∈0,T,where the geometrical property l1>0 stands for the size of Ω1 along the z-axis (see Figure 2). That is, the (chemical) equilibrium state is somewhat dependent upon the geometrical properties of the slug, as found in the paper of Bico and Quéré [1]. Thus, the expression which would match these conditions is C∞l=l2/l1+l2, since Cl=i-1 in Ωi (with i=1,2) at time t=0. Note that C∞l=1/2 is also an admissible value but only for a particular case so to speak when l1=l2, for it neglects the geometrical property of the slug, particularly, when l2→0, ∞.In regard to these arguments, one then concludes that the lesser the lengthl2>0 the higher the concentration C∞l∈0,1. It now remains to settle proper forms for γi,∞>0 (with i=1,2). First, we note that(9)limt→∞⁡l1γ1+l2γ2l1+l2=l1γ1,∞+l2γ2,∞l1+l2=l1γ1,0+l2γ2,0l1+l2on⋃i=1,2Γi;γ2-i,∞≡γ2-i,0whenl1+i=0i=0,1.Hence, these suggest to consider abarycentric-based model, namely, of the form (10)γi,∞=l1γ1,0+l2γ2,0l1+l2,onΓii=1,2,because γ1,∞=γ2,∞ (equilibrium surface tension) when t→∞. Clearly, there is no loss of generality in doing so, for not only does (10) agree well but it is also a physically acceptable expression. We are now well prepared to set forth the governing equations; this is the purpose of the next section.


## 2.4. Two-Phase Flow Phase-Field Model

In Sections2.2 and 2.3, we covered much ground on diffusive and interfacial effects. The present section takes those effects into account and models the slug self-propulsion mechanism using thephase-field approach. In a sense, we have beforehand the motion of a two-miscible-liquid system partly in contact with an immiscible one, the vapour phase region. Therefore, the starting point is naturally based on the classical Navier-Stokes equations (NS-model, say) for incompressible flow with interfacial forces coupled with the so-calledCahn-Hilliard evolution equation (CH-model, say). In the absence of gravitational field, those equations in vector invariant forms read as follows: (11)∂tρu+divρu⊗u-divTp,u-G∇ϕ=0in⋃i=1,2,vΩi×0,T,∂tϕ+divϕu-divϖ∇G=0in⋃i=1,2,vΩi×0,T,divu=0in⋃i=1,2,vΩi×0,T,where the stress tensor Tp,u, constructed with the pressure and velocity fields pr,t, ur,t reads as follows:(12)Tp,u=-pI+2μSu,whereSu=12∇u+∇uT.It is understood that Ωv=Ωv-∪Ωv+ (=Ωv±). In (11) and (12), the property constants ρ,μ>0 characterise themixed density and dynamic viscosity of the slug; both are defined in terms of the material properties of the two phases, and the smooth (phase-field) function ϕ12r,t∈0,1 is β=βv+β2-βvϕ12 (with β=ρ, μ). The function Gϕ defines the so-calledbulk chemical potential as(13)Gϕ=-ζε2divε2∇ϕ-ϕ2-1ϕ,whereϕr,t∈-1,1.In (13), ε>0[m] is acapillary width that scales with the thickness of the boundaries Γ1, Γ2; ϖ>0 [m3·s·kg−1], themobility term which determines the time scale of the Cahn-Hilliard diffusion; and ζ>0[N], themixing energy density. For more details see, for instance, the papers of Yue et al. [39] and Jacqmin [40] as well. We draw attention to the fact that ϕ12≃1+ϕ/2, and, therefore(14)ϕ12=min⁡max⁡1+ϕ2,0,1,whereϕr,t∈-1,1.In (14), min· and max· are operators which give the minimum and maximum of their respective arguments. Throughout the following discussion, the phase-field function ϕr,t∈-1,1 is used instead of ϕ12r,t∈0,1. It now remains to supplement the appropriate boundary conditions on Γ1,Γ2,Γs. At theinlet andoutlet boundaries azero pressure is imposed;that is, the system containing the domains Ω1, Ω2 is embedded in astill environment Ωv.Wetted wall boundary conditions are imposed along Γs; the contact angles are such that θ1≡θ2 (≡π/3 [rad]). The choice θk≡π/3 [rad] (with k=1,2) is justified from the fact that the advancing and receding boundaries of the miscible bidroplet system Ω1∪Ω2 are convex curves, following the works of Bico and Quéré [1, 5]. More features based partly on the field ϕr,t∈-1,1 are discussed in Section 3.


## 3. Numerical Implementation

This section furnishes the mathematical model so far discussed with numerical features to simulate the motion of the slug. Figure3 shows a simplified computational box model with the corresponding zero level set function ϕvr,t∈0,1 delineating the regions of interest.Figure 3
Sketch diagram of the computational box model used in this study paralleled with the correspondingzero level set function ϕvz,t∈0,1 to delineate the regions of interest.Recall that we are concerned with the motion of a slug so in order to adapt the solution ofClr,t with its motion one shouldremap certain property variables attached to Ω1∪Ω2 on aphase-field basis to track the evolution of this function while following the region Ω1∪Ω2. In other words, this suggests to remap the scalar fields Clr,t, φk as follows: Cl↦ϕCl=Cϕl, φk↦ϕφk=φϕ,k (with k=1,2,s). Moreover, since ϕ∈-1,1, which is equivalent to say that min⁡ϕ=-max⁡ϕ=1, the remapping ϕ↦1+ϕ/2 is first taken into account so as to obtain ϕ∈0,1, as can be easily checked by substitution. Clearly, from these remappings one observes that the convection-diffusion equation vanishes only outside the region Ω1∪Ω2. Consequently, these now result in a model which spans over ⋃i=1,2,v‍Ωi, while retaining its characteristics only in Ω1∪Ω2. Hence, invoking these mapping properties into (5) yields (15)ϕ∂tCl+divΞϕ=Fϕin⋃i=1,2,vΩi×0,T,Cϕlr,t=i-1inΩi×t=0i=1,2,-ni·Cϕlu-D21∇Cϕl=φϕ,i=-ni·CϕluonΓi×0,Ti=1,2,s,where the fields Ξϕr,t, Fϕr,t are for shortness defined by(16)Ξϕr,t=u-D21∇Cϕlin⋃i=1,2,vΩi×0,T;Fϕr,t=-∂tϕClin⋃i=1,2,vΩi×0,T.Following (15), the problem then becomes one of finding the solution of aGeneral Form PDE. Alternatively, one can formulate the convection-diffusion equation describing the scalar field Clr,t as follows:(17)∂tCl+divϕuCl-D21ϕ∇Cl=Rϕr,tin⋃i=1,2,vΩi×0,T,where the reaction term Rϕr,t appearing in the above transport reads as(18)Rϕr,t=1-ϕ∂tCl-∂tϕClin⋃i=1,2,vΩi×0,T.This is the problem of finding the solution of aTransport PDE. In COMSOL, the latter is calledTransport of Dilute Species Interface. It may be observed in passing the fact that whenever r=r,zT is a point of Ω1∪Ω2, (1-ϕ)=0 implies that ∂tϕ=0 and Cϕl=Cl, thence the original convection-diffusion equation. Contrarily, when r=r,zT is not a point of Ω1∪Ω2, (1-ϕ)=1 implies that ∂tϕ=0 and Cϕl=0, thence the restriction of the property Cl(r,t) within the domain Ω1∪Ω2. These latter two statements confirm very clearly the validity of the above underlying mapping procedures.Now, some techniques which have rendered efficient aid towards implementing this model with COMSOL are in order. Henceforth, we presume the following COMSOL interfaces are activated:the phase-field andGeneral Form PDE interfaces, respectively: the former for theNS- andCH-models and the latter for theCD-model. Decidedly, it is observed that the conditions imposed on Clr,t imply that(19)CϕlΓ2≳CϕlΓ1≥0∀t∈0,T;thence, it follows from (19), using the min· and max· operators, that(20)CϕlΓk≡max⁡Cϕlr,tfort∈0,Tifk=2;min⁡Cϕlr,tfort∈0,Tifk=1.With these features alongside, we first utilise the COMSOL min· and max· operators to construct the interfacial forces γi(CΓil) at their respective boundaries Γ1,Γ2. Then, these and other fluid properties are entered into thefluid propertiesi (with i=1,2) features. Finally, Γ1, Γ2 are selected asinitial interfaces. In regard to theGeneral Form PDE interface, (15), its coefficients and flux components, namely, ∂λΞr,t (with λ=r,z), are entered into the appropriate default nodes, the concentration in Ω2 is initialized to unity, and azero flux boundary condition is assumed from the surroundings. This being so, several simulations were then run with an array of data extracted from the papers of Bico and Quéré [1, 5]. Extracted data from those authors are tabulated in Table 1. These simulations are discussed through several cases in the following sections.Table 1
The set of property constants used to run the numerical experiments: data used for case studies 1 to 4.
Data used for case studies R (l1, l2) (ρ1, ρ2) (μ1, μ2) (γ1,0, γ2,0) [mm] [mm] [g⋅cm−3] [kg⋅m−1⋅s−1] [mN⋅m−1] 4 × 1 0 - 1 (15, 15) (1.1132, 0.95) (1.6, 1.7) × 10−2 (47.7, 20.3)


## 4. Results and Discussion

In this section we analyse the results of our simulations based on the data given in Table1. The data appearing in Table 1 are those extracted from the papers of Bico and Quéré [1, 5]. This made the need of a condensed notation imperative. Define in compact form those property constants, namely, R,ls,D21,μ>0, upon which the numerical experiments were carried out as π=defR,ls,D21,μT, and henceforward view the surface tension difference between both menisci as a function of the set of parameters π and of time t>0, respectively; that is,(21)γt,π=γ1t,π-γ2t,πfor  givenπ,t∈0,T.Let theposition andspeed/velocity,that is, Xst,π and Ust,π, respectively, of the slug be functions of the forms(22)Ξs=Ξst,γforΞs=Xs,Us;givenπ,t∈0,T,so that the corresponding differentials with respect to time write as (23)dΞsdtr,t=∂t+∂tγddγΞsr,t=LγΞsr,tforΞs=Xs,Us.From now on, we can then characterise ∂tΞs as theinstantaneous rate of change but dΞs/dγ, thenet change in the propertyΞst,γwith respect to interfacial variation and ∂tγt,π, thenet rate of change in interfacial forces with respect to time. In particular, by LγΞsr,t is meant the time derivative of Ξst,γ following the dynamics of γt,π, where Ξs=Xst,π, Ust,π. We will come back to those terms in later sections. In other words, these bring forward the key factors that do have significant effects on the evolution of the scalar fields Xst,π, Ust,π. Suppose the following baseline properties are granted: π0=(R0,ls0,D210,μ0)T, and important questions related to those properties might then be addressed. For instance, it would be fruitful to extend theoretical understanding based on the following inequalities: (24)β≤β0,β>β0forβ=R,ls,D21,μ>0.Henceforth, we primarily aim at answering two sets of questions. In the first place, we aim to answer the following questions.Question 2 (givenπ0).
What are the kinematics of the slug in steady and unsteady regimes?Question 3 (givenπ0).
How does the slug evolve whenπ0 is perturbed to π0+δπ?And, in the second place, we aim to answer the following question.Question 4 (givenπ0 and t∞,Xs,∞, resp.).
Are the positions of the slug predictable for allt∈t0,t∞?In the above last question, the scalar fieldXs,∞>0 denotes the position of the miscible fluids system after self-propulsion and t∞≤T the corresponding time. Thus, by the pairs of scalars t0,Xs,0, t∞,Xs,∞ are meant the ending points of its path in the t,Xs-plane. In the following sections we shed lights on the above three questions by advancing a series of graphical results and mechanical explanations.


## 5. The Phase-Field Function

The proposed methodological approach rests on the phase-field technique. Knowing the subsequent positions of the phase-field functionϕr,t when time evolves, one can take advantage if it calculates the subsequent position of the slug. Thus, the discussion commences by computing the scalar field ϕr,t. To do so, we use the set of property constants tabulated in Table 1.The time evolution of the phase functionϕr,t is illustrated at the sequence of times tkk=1,2,3=(5k-1)k=1,2,3 [s] in Figure 4. We draw attention that ϕr,t is graphed by the remapping ϕ↦1+ϕ/2 in order to have ϕ∈0,1 instead of ϕ∈-1,1.Figure 4
Plots of the phase-field functionϕr,t along the axis of the capillary tube at the sequence of times tkk=1,2,3=(5k-1)k=1,2,3 [s].


## 6. Case 1: Motion of the Slug versus Time

The following two cases concern the rate of penetration of the slug into the capillary tube in steady and unsteady regimes. The steady case/regime refers to the case for which the position-time graph of the slug is linear and the unsteady case/regime to other cases.


### 6.1. Case 1.1: Steady Regime

Generally, on aPoiseuille flow hypothesis together with the assertion that the slug issufficiently long, ls≫1 say, its path results in an integral expression of the form(25)Xst,γ=∫γ1-γ2R4μ1l1+μ2l2dt+Const.for  givenπ,t∈0,T.To settle (25), we have taken into account the following statements [1]: the Poiseuille flow induced within the region Ω1∪Ω2 is of the form ur=δpR2-r2/4μsls, where the Laplacian pressure exerted at Γ1relative to that exerted at Γ2 reads as δp=2γ1-γ2/R. The velocity field ur is averaged in Ω1∪Ω2, and the resulting expression equated to dXs/dt, yielding dXs/dt=Rγ1-γ2/4μ1l1+μ2l2, since μls=μ1l1+μ2l2, thence the antiderivative of (25).Consequently, it follows that there exists a functionalχt,π such that Xst,γ=∫χt,πdt+Const. for given π=R,ls,D21,μT, t>0. Numerical simulation was carried out for the steady case in the hope of establishing the character of χt,π. Theearly-time results confirm that χ=χπ, and the position of the slug obeys the following linear law at early-time:(26)Xst,γ≈χπ·t-t0,χπ≈Const.  fort∈t0,t10⊂0,T,where t10=10 [s], as can be easily seen in Figure 5. To graph Figure 5, the COMSOL Cut Line data set is used to create a line parallel to the axis of the capillary tube, having its ending points at r,z=R/2,0, R/2,L, where L≫ls denotes the length of the capillary tube along the z-axis. This allows to track the subsequent positions occupied by the advancing and receding boundaries,that is, Γ1, Γ2, since it will cut at every instant of time t∈0,T these boundaries at precisely two and only two (boundary) points r=R/2,zΓktT∈Γk (with k=1,2). Setting Xs=zΓ1t, we then construct a time-position sequence tk,Xstk,πk∈0,10=k,Xs,kk∈0,10. Representing these data points by circular markers (∘), the corresponding approximation by solid line, and that of Bico and Quéré [5] by dotted line, one obtains Figure 5. Henceforth, whenever the study is carried out with respect to time, it is to be understood that the plotting referred to in the t,Xs-plane is based on that approach, unless the study is carried out parametrically, in which case we shall then draw attention upon.Figure 5
Position of the slug as function of time,Xst≈1.5t-0.15 [mm]. The radius of the capillary tube is 0.4 [mm], and the successive lengths are l1=16 [mm] and l2=11 [mm]. Note that the experimental result gives Us≈1.55 [mm⋅s−1], and ours reads as Us≡∂tXs≈1.50 [mm⋅s−1]. The numerical result (solid line) is compared with that of Bico and Quéré [5] (dotted line).In the sequel, we will show through numerical experiments thatχπ≈Rγ1-γ2/4μ1l1+μ2l2. Thus, by (26) it would appear that after some time t>t0, the slug will have travelled at aconstant velocity Us≡χπ a distance Xst,γ≈χπ·t-t0. Clearly, (26) appears to show that Xs→∞ whenever t→∞, a situation which is reasonable mathematically but unreasonable physically; this is simply because the mathematical description using the linear approximation Xst,γ≈χπ·t-t0 will at some time t>0 (unsteady case) fail to give physical agreement. In that case a nonlinear law (as will be shown in the sequel to exist) should be sought.The above statements suggest that there exists a scalarXs∗∈0,∞, characterising the distance at which the slug stops, such that limt→∞Xst,γ→Xs∗, since we know from the nature of the problem that the difference in interfacial forces, which was initially very high,that is, γ1,0-γ2,0 (≡γ0), owing to the initial conditions Clr,t=0=i-1 in Ωi (with i=1,2), will have decreased to a vanishingly small value when the equilibrium of Clr,t has been established. Thus, this demonstrates the existence of a temporal (π-dependant) property t∞∈0,T such that mechanical equilibrium is attained via interfacial effects;namely,(27)γ1t,π≈γ2t,π,whenevert≳t∞∈0,T;πgiven.Clearly, there is no loss of generality in assigning, for instance, this particular case to an initially prescribed set, say β≡β0, where β=R,ls,D21,μ. Based on the said facts, we ran the numerical simulation with a large final time T=100 [s] to demonstrate that (27) is satisfied.Surprisingly, we observed with a slug bearing the characteristicsls=30 [mm], μk≈17 [mPa·s] (with k=1,2), and D21≈10-9 [m2·s]—the diffusion constant for alcohol fluid 1 in water fluid 2 [41, 42]—that(28)γt∞,π≲5%ofmin⁡γ1,0,γ2,0for  givenπ,t∞≳100s.The graphical evolution of the surface tension function γt,π with respect to time t∈t0,t∞ is shown in Figure 6. When the surface tensions γit;π (i=1,2) are plotted on the same t,γ-plane, one obtains Figure 7, where curved arrows appearing in that plane denote the direction of convergence of these functions with respect to time.Figure 6
Surface tension difference versus time:γt,π [N⋅m−1]. The radius of the capillary tube is 0.4 [mm]; the length of the slug is such that ls=l1+l2≈15ls+15ls [mm] (with both the scalar fields ls, π keptfixed). Observe that γt∞,π≲5% of minγ1,0,γ2,0 for t∞≳100 [s].Figure 7
The time evolution of the surface tensionsγit;π (i=1, 2).A very interesting feature is observed from Figure7. We have γ1t;π≈γ2t;π≈3.2×10-2 [N·m−1] at time t≈100 [s]. Substituting the values tabulated in Table 1 into the barycentric-model equation (10) for the pairs of scalar fields l1,l2, γ1,0,γ2,0, one arrives at γi,∞≈3.3×10-2 [N·m−1] (with i=1,2), thereby validating the form of (10).Out of the graphical results displayed in Figures6 and 7, it would appear that (26) does not hold in both the early- and late-time regimes, more precisely, in thesteady (orearly-time) regime as well as in theunsteady (orlate-time) regime. The next section discusses the dynamics of the slug in the unsteady regime.


### 6.2. Case 1.2: Unsteady Regime

When the period over which the numerical experiments were extended tot0,≲150 [s], the subsequent positions of the slug graphed exponentially, as can be easily shown in Figure 8. Hence, the path of the slug could fairly well be described by anexponential-type equation of the form(29)Xst,π≈η11-exp⁡-η2t-t0,for  givenπ,t∈0,T.In (29), the fields η1, η2>0 are π-dependent scalars;that is, η1=η1π and η2=η2π. Furthermore, since the slug movessteadily only for a short period of time 0,t10 (roughly, t10≈ 10 [s]), consequently, there does necessarily exist Xs,∞, Us>0 such that, for given π=R,ls,D21,μT,(30)limt→t∞⁡Xst,π⟶Xs,∞,limt→t10⁡∂tXst,π⟶Us,limt→t∞⁡∂tXst,π⟶0,where 0≤t0<t10<t∞≤T. Here, it is to be understood that Xs,∞>0 is theplace at which the slugstops whereas Us>0 simply characterises itssteady regime. Also, t0,t10,t∞≥0 stand, respectively, for the initial time t0=0 [s], the time t10=10 [s], and the time t∞≤T [s] at which the slug stops. But T>0 designates the period over which the numerical experiments are carried out. With these characteristics at hand, the path of the slug then writes as(31)Xst,π≈Xs,∞1-exp⁡-UsXs,∞t-t0≡Ust-t0-∑k=2∞UskXs,∞k-1k!t0-tk≡Xs,0t,π+xst,πfor  givenπ,t∈0,T.From this, one concludes that the position Xst,π of the slug is merely a superposition arising from two contributions: thesteady andunsteady regimes; the former is characterised by Xs,0t,π and the latter by xst,π, respectively. Following Figure 5, it is obvious that the equality Xs,0=Xst,π holds only for t∈t0,10 [s]. Thus, for early times, the position of the slug can be approximated as(32)Xst,π≈Ust-t0for  givenπ,t∈0,t10⊂0,T,which is precisely the same expression in the steady case, Us>0 being put for η1·η2>0. It now remains to demonstrate that Us (≈η1·η2) coincides with the result of Bico and Quéré [5].Figure 8
Slug position versus time, described by (31) and compared with numerical results (empty squares (□) for ls=30 [mm] and empty circles (∘) for ls=60 [mm]).For two typical slugs of lengthsls=30, 60 [mm], respectively, we combined (29) with numerical data in the least-squares sense. We found η1,η2≈50.0542,0.0306 for ls=30 [mm] and η1,η2≈49.8052,0.0129 for ls=60 [mm]. Thence, with the former, we deduce Us≈1.5312 [mm·s-1] and with the latter Us≈0.6425 [mm·s−1]. All these quantities are arranged carefully in Table 2. Figure 8 is intended to exhibit these features graphically using (31).Table 2
Numerical values for (η1, η2), data for (29).
Results of case study 1.2 (η1, η2) l s = 30 [mm] l s = 60 [mm] η 1 [mm] 50.0542 49.8052 η 2 [s-1] 0.0306 0.0129 U s ≈ η 1 ⋅ η 2 [mm⋅s-1] 1.5312 0.6425No doubt, it is no exaggeration to say that our model, verified numerically in those two cases, namely, forls=30, 60 [mm], as depicted in Figure 8, may be extended to analyse the movement of small volumes of viscous liquids in capillary tubes of square, rectangular cross sections and so forth. Hence, given the set Xs,∞,Us, that pair of scalar fields designating the point at which the slug stops for an initially given speed, one can predicta priori the position as well as the rate of penetration of the slug into the capillary tube. Further, if the focus is onrapid motion, one should evidently consider a small slug, instead of a big one. Contrarily, if the primary goal is to achieve large penetrating length, one might endeavour to take a large slug; the trueness of these two statements is reflected in Figure 8. The following sections enable further physical relationships to be grasped.


## 6.1. Case 1.1: Steady Regime

Generally, on aPoiseuille flow hypothesis together with the assertion that the slug issufficiently long, ls≫1 say, its path results in an integral expression of the form(25)Xst,γ=∫γ1-γ2R4μ1l1+μ2l2dt+Const.for  givenπ,t∈0,T.To settle (25), we have taken into account the following statements [1]: the Poiseuille flow induced within the region Ω1∪Ω2 is of the form ur=δpR2-r2/4μsls, where the Laplacian pressure exerted at Γ1relative to that exerted at Γ2 reads as δp=2γ1-γ2/R. The velocity field ur is averaged in Ω1∪Ω2, and the resulting expression equated to dXs/dt, yielding dXs/dt=Rγ1-γ2/4μ1l1+μ2l2, since μls=μ1l1+μ2l2, thence the antiderivative of (25).Consequently, it follows that there exists a functionalχt,π such that Xst,γ=∫χt,πdt+Const. for given π=R,ls,D21,μT, t>0. Numerical simulation was carried out for the steady case in the hope of establishing the character of χt,π. Theearly-time results confirm that χ=χπ, and the position of the slug obeys the following linear law at early-time:(26)Xst,γ≈χπ·t-t0,χπ≈Const.  fort∈t0,t10⊂0,T,where t10=10 [s], as can be easily seen in Figure 5. To graph Figure 5, the COMSOL Cut Line data set is used to create a line parallel to the axis of the capillary tube, having its ending points at r,z=R/2,0, R/2,L, where L≫ls denotes the length of the capillary tube along the z-axis. This allows to track the subsequent positions occupied by the advancing and receding boundaries,that is, Γ1, Γ2, since it will cut at every instant of time t∈0,T these boundaries at precisely two and only two (boundary) points r=R/2,zΓktT∈Γk (with k=1,2). Setting Xs=zΓ1t, we then construct a time-position sequence tk,Xstk,πk∈0,10=k,Xs,kk∈0,10. Representing these data points by circular markers (∘), the corresponding approximation by solid line, and that of Bico and Quéré [5] by dotted line, one obtains Figure 5. Henceforth, whenever the study is carried out with respect to time, it is to be understood that the plotting referred to in the t,Xs-plane is based on that approach, unless the study is carried out parametrically, in which case we shall then draw attention upon.Figure 5
Position of the slug as function of time,Xst≈1.5t-0.15 [mm]. The radius of the capillary tube is 0.4 [mm], and the successive lengths are l1=16 [mm] and l2=11 [mm]. Note that the experimental result gives Us≈1.55 [mm⋅s−1], and ours reads as Us≡∂tXs≈1.50 [mm⋅s−1]. The numerical result (solid line) is compared with that of Bico and Quéré [5] (dotted line).In the sequel, we will show through numerical experiments thatχπ≈Rγ1-γ2/4μ1l1+μ2l2. Thus, by (26) it would appear that after some time t>t0, the slug will have travelled at aconstant velocity Us≡χπ a distance Xst,γ≈χπ·t-t0. Clearly, (26) appears to show that Xs→∞ whenever t→∞, a situation which is reasonable mathematically but unreasonable physically; this is simply because the mathematical description using the linear approximation Xst,γ≈χπ·t-t0 will at some time t>0 (unsteady case) fail to give physical agreement. In that case a nonlinear law (as will be shown in the sequel to exist) should be sought.The above statements suggest that there exists a scalarXs∗∈0,∞, characterising the distance at which the slug stops, such that limt→∞Xst,γ→Xs∗, since we know from the nature of the problem that the difference in interfacial forces, which was initially very high,that is, γ1,0-γ2,0 (≡γ0), owing to the initial conditions Clr,t=0=i-1 in Ωi (with i=1,2), will have decreased to a vanishingly small value when the equilibrium of Clr,t has been established. Thus, this demonstrates the existence of a temporal (π-dependant) property t∞∈0,T such that mechanical equilibrium is attained via interfacial effects;namely,(27)γ1t,π≈γ2t,π,whenevert≳t∞∈0,T;πgiven.Clearly, there is no loss of generality in assigning, for instance, this particular case to an initially prescribed set, say β≡β0, where β=R,ls,D21,μ. Based on the said facts, we ran the numerical simulation with a large final time T=100 [s] to demonstrate that (27) is satisfied.Surprisingly, we observed with a slug bearing the characteristicsls=30 [mm], μk≈17 [mPa·s] (with k=1,2), and D21≈10-9 [m2·s]—the diffusion constant for alcohol fluid 1 in water fluid 2 [41, 42]—that(28)γt∞,π≲5%ofmin⁡γ1,0,γ2,0for  givenπ,t∞≳100s.The graphical evolution of the surface tension function γt,π with respect to time t∈t0,t∞ is shown in Figure 6. When the surface tensions γit;π (i=1,2) are plotted on the same t,γ-plane, one obtains Figure 7, where curved arrows appearing in that plane denote the direction of convergence of these functions with respect to time.Figure 6
Surface tension difference versus time:γt,π [N⋅m−1]. The radius of the capillary tube is 0.4 [mm]; the length of the slug is such that ls=l1+l2≈15ls+15ls [mm] (with both the scalar fields ls, π keptfixed). Observe that γt∞,π≲5% of minγ1,0,γ2,0 for t∞≳100 [s].Figure 7
The time evolution of the surface tensionsγit;π (i=1, 2).A very interesting feature is observed from Figure7. We have γ1t;π≈γ2t;π≈3.2×10-2 [N·m−1] at time t≈100 [s]. Substituting the values tabulated in Table 1 into the barycentric-model equation (10) for the pairs of scalar fields l1,l2, γ1,0,γ2,0, one arrives at γi,∞≈3.3×10-2 [N·m−1] (with i=1,2), thereby validating the form of (10).Out of the graphical results displayed in Figures6 and 7, it would appear that (26) does not hold in both the early- and late-time regimes, more precisely, in thesteady (orearly-time) regime as well as in theunsteady (orlate-time) regime. The next section discusses the dynamics of the slug in the unsteady regime.


## 6.2. Case 1.2: Unsteady Regime

When the period over which the numerical experiments were extended tot0,≲150 [s], the subsequent positions of the slug graphed exponentially, as can be easily shown in Figure 8. Hence, the path of the slug could fairly well be described by anexponential-type equation of the form(29)Xst,π≈η11-exp⁡-η2t-t0,for  givenπ,t∈0,T.In (29), the fields η1, η2>0 are π-dependent scalars;that is, η1=η1π and η2=η2π. Furthermore, since the slug movessteadily only for a short period of time 0,t10 (roughly, t10≈ 10 [s]), consequently, there does necessarily exist Xs,∞, Us>0 such that, for given π=R,ls,D21,μT,(30)limt→t∞⁡Xst,π⟶Xs,∞,limt→t10⁡∂tXst,π⟶Us,limt→t∞⁡∂tXst,π⟶0,where 0≤t0<t10<t∞≤T. Here, it is to be understood that Xs,∞>0 is theplace at which the slugstops whereas Us>0 simply characterises itssteady regime. Also, t0,t10,t∞≥0 stand, respectively, for the initial time t0=0 [s], the time t10=10 [s], and the time t∞≤T [s] at which the slug stops. But T>0 designates the period over which the numerical experiments are carried out. With these characteristics at hand, the path of the slug then writes as(31)Xst,π≈Xs,∞1-exp⁡-UsXs,∞t-t0≡Ust-t0-∑k=2∞UskXs,∞k-1k!t0-tk≡Xs,0t,π+xst,πfor  givenπ,t∈0,T.From this, one concludes that the position Xst,π of the slug is merely a superposition arising from two contributions: thesteady andunsteady regimes; the former is characterised by Xs,0t,π and the latter by xst,π, respectively. Following Figure 5, it is obvious that the equality Xs,0=Xst,π holds only for t∈t0,10 [s]. Thus, for early times, the position of the slug can be approximated as(32)Xst,π≈Ust-t0for  givenπ,t∈0,t10⊂0,T,which is precisely the same expression in the steady case, Us>0 being put for η1·η2>0. It now remains to demonstrate that Us (≈η1·η2) coincides with the result of Bico and Quéré [5].Figure 8
Slug position versus time, described by (31) and compared with numerical results (empty squares (□) for ls=30 [mm] and empty circles (∘) for ls=60 [mm]).For two typical slugs of lengthsls=30, 60 [mm], respectively, we combined (29) with numerical data in the least-squares sense. We found η1,η2≈50.0542,0.0306 for ls=30 [mm] and η1,η2≈49.8052,0.0129 for ls=60 [mm]. Thence, with the former, we deduce Us≈1.5312 [mm·s-1] and with the latter Us≈0.6425 [mm·s−1]. All these quantities are arranged carefully in Table 2. Figure 8 is intended to exhibit these features graphically using (31).Table 2
Numerical values for (η1, η2), data for (29).
Results of case study 1.2 (η1, η2) l s = 30 [mm] l s = 60 [mm] η 1 [mm] 50.0542 49.8052 η 2 [s-1] 0.0306 0.0129 U s ≈ η 1 ⋅ η 2 [mm⋅s-1] 1.5312 0.6425No doubt, it is no exaggeration to say that our model, verified numerically in those two cases, namely, forls=30, 60 [mm], as depicted in Figure 8, may be extended to analyse the movement of small volumes of viscous liquids in capillary tubes of square, rectangular cross sections and so forth. Hence, given the set Xs,∞,Us, that pair of scalar fields designating the point at which the slug stops for an initially given speed, one can predicta priori the position as well as the rate of penetration of the slug into the capillary tube. Further, if the focus is onrapid motion, one should evidently consider a small slug, instead of a big one. Contrarily, if the primary goal is to achieve large penetrating length, one might endeavour to take a large slug; the trueness of these two statements is reflected in Figure 8. The following sections enable further physical relationships to be grasped.


## 7. Case 2: Motion of the Slug versus Geometrical and Physical Properties

The position of the slug is such thatXs=Xst,π. So far, we studied this position function with respect to time, keeping the elements of the parameter vector π=(R,ls,D21,μ)T fixed at those values which are tabulated in Table 1. Here it is that this function is studied parametrically. We agree that parametric studies are best undertaken when the property in question is time-invariant. We have, on the one hand, limt→t∞Xst,π→Xs,∞=Const. and limt→t10∂tXst,π→Us=Const.; see (30). And, on the other hand, the slug moves according to Poiseuille flow only in the early-time regime. As a consequence, it seems, therefore, legitimate to carry out the parametric studies in the early-time regime, which are what we proposed to consider hereinafter.Letμs>0 designate the dynamic viscosity of the slug (mixture of the two miscible liquid droplets). On abarycentric-based basis we may then consider μ·ls=μ1l1+μ2l2, in which case the velocity Ust,π of the slug, following (29), writes as(33)Ust,π≈Rγ1-γ24μ1l1+μ2l2≈Rγ1-γ24μlsfor  givenπ,t∈0,T.We understand that (33) is the derivative of (25) with respect to t∈0,T;that is, ∂tXs=Ust,π. If (33) holds true, so are the following proportionalities: Ust,π∝R, μ-1, and ls-1>0. These are what we are going to show through a series of numerical experiments in the following sections.


### 7.1. Case 2.1: Motion of the Slug versus the Radius of the Capillary Tube

In this section the motion of the slug is studied in a series of capillary tubes having the following radii:Rk=2k×10-1 [mm] (with k∈1,5). The properties of the slug are unaltered. We aim at establishing an algebraic expression relating the velocity Us of the slug to the radius R of the capillary tube. For every k∈1,5, a numerical experiment is carried out and the velocity Usk>0 of the slug is calculated. When, for each k∈1,5, the ratio Usk/Rk is analysed, it was found that(34)UskRk≈2.4×10-2±5%s-1∀k∈1,5,from which we conclude that Us∝R. The corresponding pairs of values (Rk,Usk)k=15 are plotted on a R,Us-plane; see Figure 9. Consequently, the dependence of the functional χπ finally restricted to ls,D21,μ≡Consts evolves linearly with respect to the radius R of the capillary tube according to the following linear relationship:(35)χπ≈KRR+KR,0,whereKR,KR,0=Consts.It is understood that the scalar fields KR, KR,0 are partially π-dependant constants (see the character of (33) to understand this point). In other words, an increase in R by δR>0 results in a linear increase in Us>0, conversely, which manifests the nature of the size of the capillary tube. The next section examines the influence of the geometrical property ls>0 on the kinematics of the slug.Figure 9
Plots of the sequence of points(Rk,Usk)k=15. The linear expression describing the solid line is Us≈2.4R-0.08. The properties of the slug are l1=16 [mm], l2=16.5 [mm], and μk=17 [mPa⋅s] (with k=1, 2). Other properties used are those given in Table 1.


### 7.2. Case 2.2: Motion of the Slug versus Its Length

We herein explore the effect of the geometrical propertyls>0 (size of the slug) on the slug motion. We ran for the sequence ls=30,40,…,100 [mm] (with lk=ls/2, k=1,2) several numerical experiments, keeping R,D21,μ≡Consts. The accompanying table (Table 3) compares the typical parameters for the slug self-propulsion and its laboratory analogs extracted from the paper of Bico and Quéré [5].Table 3
Comparison ofXst,π=χπ⋅t-t0 in the early-time regime.
Results of case study 2.2 l s X s t U s ≡ χ π Bico  and Quéré [5] [mm] [mm] [mm⋅s-1] [mm⋅s-1] 30 1.7 t - 0.29 1.7 1.6 40 1.1 t + 0.76 1.1 1.0 50 0.89 t + 0.76 0.89 0.8 60 0.77 t + 0.82 0.77 0.7 70 0.69 t - 0.82 0.69 0.65 80 0.57 t + 0.88 0.57 0.6 90 0.53 t - 0.79 0.53 0.5 100 0.48 t + 0.8 0.48 0.4Figure10 gives some graphical illustrations on this basis. On a careful perusal of Table 3 and Figure 10, respectively, one concludes that the agreement between simulations and experiments is good. As a result, the present study puts forward an important point: an increase in ls>0 (by δls=10 [mm], say) weakens the surface tension difference, that is, γt,π. Hence, one concludes that a slug of length ls movesfaster than the one whose length slightly differs, for instance, by δls>0. Furthermore, perusing Figure 10, we conclude that the dependence of the functional χπ restricted to R,D21,μ=Consts on ls>0 evolves on a log-log scale according to(36)log⁡χπ≈-Klslog⁡ls+Kls,0,whereKls,Kls,0=Consts.As pointed out so far, it is understood that the scalar fields Kls, Kls,0 are partially π-dependant constants. From Figures 10 and 11, it may be remarked in passing that(37)∂tXsls=30≥∂tXsls≥30for  givenπ,t∈0,T.Converting every log-log size-speed (log⁡ls,log⁡Us) pair into the usual (ls,Us) pair and, further, letting Usnum, Usexp>0 be the corresponding numerical, experimental values, the tabulation of all these property constants yields Table 4.Table 4
Numerical and experimental magnitudes of the speedsUsnum, Usexp>0, respectively, of the slug versus the size ls>0 of the slug.
Results of case study 2.2 l s  [mm] 30 40 50 60 70 80 90 100 U s n u m  [mm⋅s-1] 1.7 1.1 0.89 0.77 0.69 0.57 0.53 0.48 U s e x p   [mm⋅s-1] 1.5 1.0 0.8 0.7 0.65 0.6 0.5 0.4Figure 10
Position of the slug restricted tols=30,50,70,90 [mm] (with lk=ls/2,  k=1,2) as function of time (and of the property π): Xst,π. The radius of the capillary tube is 0.4 [mm].Figure 11
Slug velocity versus its lengthls [mm], logUs≈-1.0061log⁡ls+3.8724. The radius of the capillary tube is 0.4 [mm], and the successive lengths are l1=16 [mm] and l2=16.5 [mm]. The numerical result (solid line) is compared with that of Bico and Quéré [5] (dotted line). Note that from the graphical result of Bico and Quéré [5] we deduce log⁡Us≈-0.97297log⁡ls+3.652.


### 7.3. Case 2.3: Motion of the Slug versus Diffusivity

This section discusses the effects of diffusivity constantD21>0. A way of interpreting the latter may be understood on the following basis. At leading order, we deduced from the constitutive equation of Clr,t the following scaling (diffusion length): lt≈D21t, which gave insight of how far the concentration has propagated inside the slug in time t>0 for a given diffusivity constant D21>0. In Figure 12 is graphed on a log-log scale the speed Us>0 of the slug versus the diffusivity constant D21>0. When every log-log diffusivity-speed (log⁡D21,log⁡Us) pair is converted into the usual (D21,Us) pair and, further, Us=Usnum>0 to stress that the dependence of the speed of the slug on the diffusivity has not been studied numerically, the tabulation of these property constants yields Table 5.Table 5
Numerical magnitudes of the speedUsnum>0 of the slug versus the diffusivity D21>0 of the slug.
Results of case study 2.3 D 21 [m2⋅s-1] 1 0 - 6 1 0 - 7 1 0 - 8 1 0 - 9 U s n u m [mm⋅s-1] 1.7 1.8 1.9 2.0Figure 12
Slug velocity versus diffusivityD21 [m2⋅s−1] in logarithmic scale; solid line is represented by log⁡Us≈-0.023522log⁡D21+0.20714. The radius of the capillary tube is 0.4 [mm], and the successive lengths are l1=16 [mm] and l2=16.5 [mm].As a result, incrementing the diffusivity constantD21>0 by δD21>0 yields(38)γt,πD21+δD21<γt,πD21for  givenπ,t∈0,T.On the whole, augmenting D21>0 therefore compels the propulsion force, which is proportional to the surface tension difference γt,π, to decrease in magnitude. As a consequence, the transport coefficient D21>0 therefore does have bearing on the motion of the slug.In other words, augmenting the magnitude of the diffusivity constantD21>0 results in not only diminishing the magnitude of the slug velocity Us>0, but also causing the propulsion force to vanish/die earlier. As shown in Figure 12, one clearly sees that the π-dependent functional χπ now restricted to R,ls,μ=Consts evolves logarithmically with the diffusion constant D21>0 according to the following log-log relationship:(39)log⁡χπ≈-KD21log⁡D21+KD21,0whereKD21,KD21,0=Consts.As pointed out so far, it is understood that the scalar fields KD21, KD21,0 are partially π-dependant constants. This is an important feature which has not been considered in any of the two papers of Bico and Quéré [1, 5]. It may be remarked further that (39) tends to prove that the logarithm of the speed of the slug log⁡Us (since χπ≡Us) is proportional to the logarithm of the diffusivity constant log⁡D21. Moreover, since lt≈D21t (that diffusion length which provides a measure of how far the chemical concentration has propagated into the slug in the z-direction by diffusion in time t>0), the existence of a reference time τref>0 can be assumed such that(40)Ordlog⁡Us~Ordlog⁡dldt~log⁡D21τref+Const.,from which the shape of (39) at once comes to light. Physically speaking, if the diffusion constant D21>0 is to furnish us with a means of strengthening the self-propulsion mechanism, it is evident that it should below in magnitude. Hence, this establishes a relation between the diffusion constant D21>0 and the speed Us>0 of the slug. The next section considers the speed Us>0 of the slug versus the viscosity constant μ>0.


### 7.4. Case 2.4: Motion of Slug versus Viscosity

The nature of the movement of the slug with the character of the viscous force is here considered by establishing a relationship relating the speedUs>0 of the slug to the viscosity constant μ>0. To effect this purpose we ran several simulations, changing the magnitude of the property μ>0 at each run while restricting other property constants to those values appearing in Table 1. Further, at each run the pair of scalar fields μ,χπ is noted. In Figures 13 and 14 are plotted the pairs of scalar fields μi,χπi=1,2 on a log-log scale. As one can see from Figures 13 and 14, the dependence of the functional χπ finally restricted to R,ls,D21≡Consts on μi>0 (with i=1, 2) evolves on a log-log scale according to(41)log⁡χπ≈-Kμlog⁡μ+Kμ,0whereKμ,Kμ,0=Consts.As pointed out so far, it is understood that the scalar fields Kμ, Kμ,0 are partially π-dependant constants. Moreover, because this expression of Us>0 in terms of D21>0 is algebraically similar to that established in the foregoing section, our thoughts are thereby led to regard a decrease in μ>0 resulting in an exponential increase in Us>0, which manifests the nature of viscous effect.Figure 13
Slug velocity versus viscosityμ2 [mPa⋅s] in logarithmic scale; numerical result, approximated by log⁡Us≈-0.63452log⁡μ2+2.6509, is graphed in solid line; experimental result extracted from earlier work approximated by logUs≈-0.57214logμ2+2.39 is in dotted line. The radius of the capillary tube is 0.51 [mm], and the successive lengths are l1=16 [mm] and l2=16.5 [mm]; also, the viscosity μ1=17 [mPa⋅s] is kept constant.Figure 14
Slug velocity versus viscosityμ1 [mPa⋅s] in logarithmic scale; numerical result, approximated by logUs≈-0.42437log⁡μ1-2.0723, is graphed in solid line. The radius of the capillary tube is 0.51 [mm], and the successive lengths are l1=16 [mm] and l2=16.5 [mm]; also, the viscosity μ2=17 [mPa⋅s] is kept constant.To the informed reader the lawUs∝μ-1 is obvious whenever the Poiseuille law is pronounced in the mathematical development but, to the uninformed reader, the law Us∝μ-1/2 is unobvious. Thus, before advancing further, it might be interesting to justify this point.Since, in the earlier section (see Section7.2), we showed that Us∝ls-1 and, from Poiseuille flow, we have Us∝μ-1, it is plain that(42)Us∝1μ1α·μ21-αifμ=μ1,μ1=μ2;α∈0,1∝1μ1αifμ2=Const.;α∈0,1∝1μ21-αifμ1=Const.;α∈0,1.Equation (42) states that the dependence of the slug speed on the dynamic viscosity μi>0 of the i-fluid (with i=1, 2) is such that Us∝μi-β for some real β∈0,1. Moreover, it also states that if such dependence is known for one fluid, the other can be deduced, since whenever Us∝μ1-β necessarily Us∝μ21-β. Otherwise the condition Us∝μ-1 will be violated. Taking advantage of the properties of logarithms to convert the second and third lines of (42) into linear relationships, one obtains(43)log⁡Us=κi·log⁡μi+Const.∀i∈1,2,where∑i=1,2κi=1.Through numerical experiments, we deduced κ1,κ2=-0.42437,-0.63452. Since ∑i=1,2‍κi=1.05889, it follows that ∑i=1,2‍κi≈1, hence satisfying (43) and, therefore, (42). We have thus justified the results appearing in Figures 13 and 14.Having studied the motion of the slug under variations of the property vectorπ=(R,ls,D21,μ)T, we propose in the next section an algebraic expression, characteristics of the nature of the movement of the slug, on this ground.


## 7.1. Case 2.1: Motion of the Slug versus the Radius of the Capillary Tube

In this section the motion of the slug is studied in a series of capillary tubes having the following radii:Rk=2k×10-1 [mm] (with k∈1,5). The properties of the slug are unaltered. We aim at establishing an algebraic expression relating the velocity Us of the slug to the radius R of the capillary tube. For every k∈1,5, a numerical experiment is carried out and the velocity Usk>0 of the slug is calculated. When, for each k∈1,5, the ratio Usk/Rk is analysed, it was found that(34)UskRk≈2.4×10-2±5%s-1∀k∈1,5,from which we conclude that Us∝R. The corresponding pairs of values (Rk,Usk)k=15 are plotted on a R,Us-plane; see Figure 9. Consequently, the dependence of the functional χπ finally restricted to ls,D21,μ≡Consts evolves linearly with respect to the radius R of the capillary tube according to the following linear relationship:(35)χπ≈KRR+KR,0,whereKR,KR,0=Consts.It is understood that the scalar fields KR, KR,0 are partially π-dependant constants (see the character of (33) to understand this point). In other words, an increase in R by δR>0 results in a linear increase in Us>0, conversely, which manifests the nature of the size of the capillary tube. The next section examines the influence of the geometrical property ls>0 on the kinematics of the slug.Figure 9
Plots of the sequence of points(Rk,Usk)k=15. The linear expression describing the solid line is Us≈2.4R-0.08. The properties of the slug are l1=16 [mm], l2=16.5 [mm], and μk=17 [mPa⋅s] (with k=1, 2). Other properties used are those given in Table 1.


## 7.2. Case 2.2: Motion of the Slug versus Its Length

We herein explore the effect of the geometrical propertyls>0 (size of the slug) on the slug motion. We ran for the sequence ls=30,40,…,100 [mm] (with lk=ls/2, k=1,2) several numerical experiments, keeping R,D21,μ≡Consts. The accompanying table (Table 3) compares the typical parameters for the slug self-propulsion and its laboratory analogs extracted from the paper of Bico and Quéré [5].Table 3
Comparison ofXst,π=χπ⋅t-t0 in the early-time regime.
Results of case study 2.2 l s X s t U s ≡ χ π Bico  and Quéré [5] [mm] [mm] [mm⋅s-1] [mm⋅s-1] 30 1.7 t - 0.29 1.7 1.6 40 1.1 t + 0.76 1.1 1.0 50 0.89 t + 0.76 0.89 0.8 60 0.77 t + 0.82 0.77 0.7 70 0.69 t - 0.82 0.69 0.65 80 0.57 t + 0.88 0.57 0.6 90 0.53 t - 0.79 0.53 0.5 100 0.48 t + 0.8 0.48 0.4Figure10 gives some graphical illustrations on this basis. On a careful perusal of Table 3 and Figure 10, respectively, one concludes that the agreement between simulations and experiments is good. As a result, the present study puts forward an important point: an increase in ls>0 (by δls=10 [mm], say) weakens the surface tension difference, that is, γt,π. Hence, one concludes that a slug of length ls movesfaster than the one whose length slightly differs, for instance, by δls>0. Furthermore, perusing Figure 10, we conclude that the dependence of the functional χπ restricted to R,D21,μ=Consts on ls>0 evolves on a log-log scale according to(36)log⁡χπ≈-Klslog⁡ls+Kls,0,whereKls,Kls,0=Consts.As pointed out so far, it is understood that the scalar fields Kls, Kls,0 are partially π-dependant constants. From Figures 10 and 11, it may be remarked in passing that(37)∂tXsls=30≥∂tXsls≥30for  givenπ,t∈0,T.Converting every log-log size-speed (log⁡ls,log⁡Us) pair into the usual (ls,Us) pair and, further, letting Usnum, Usexp>0 be the corresponding numerical, experimental values, the tabulation of all these property constants yields Table 4.Table 4
Numerical and experimental magnitudes of the speedsUsnum, Usexp>0, respectively, of the slug versus the size ls>0 of the slug.
Results of case study 2.2 l s  [mm] 30 40 50 60 70 80 90 100 U s n u m  [mm⋅s-1] 1.7 1.1 0.89 0.77 0.69 0.57 0.53 0.48 U s e x p   [mm⋅s-1] 1.5 1.0 0.8 0.7 0.65 0.6 0.5 0.4Figure 10
Position of the slug restricted tols=30,50,70,90 [mm] (with lk=ls/2,  k=1,2) as function of time (and of the property π): Xst,π. The radius of the capillary tube is 0.4 [mm].Figure 11
Slug velocity versus its lengthls [mm], logUs≈-1.0061log⁡ls+3.8724. The radius of the capillary tube is 0.4 [mm], and the successive lengths are l1=16 [mm] and l2=16.5 [mm]. The numerical result (solid line) is compared with that of Bico and Quéré [5] (dotted line). Note that from the graphical result of Bico and Quéré [5] we deduce log⁡Us≈-0.97297log⁡ls+3.652.


## 7.3. Case 2.3: Motion of the Slug versus Diffusivity

This section discusses the effects of diffusivity constantD21>0. A way of interpreting the latter may be understood on the following basis. At leading order, we deduced from the constitutive equation of Clr,t the following scaling (diffusion length): lt≈D21t, which gave insight of how far the concentration has propagated inside the slug in time t>0 for a given diffusivity constant D21>0. In Figure 12 is graphed on a log-log scale the speed Us>0 of the slug versus the diffusivity constant D21>0. When every log-log diffusivity-speed (log⁡D21,log⁡Us) pair is converted into the usual (D21,Us) pair and, further, Us=Usnum>0 to stress that the dependence of the speed of the slug on the diffusivity has not been studied numerically, the tabulation of these property constants yields Table 5.Table 5
Numerical magnitudes of the speedUsnum>0 of the slug versus the diffusivity D21>0 of the slug.
Results of case study 2.3 D 21 [m2⋅s-1] 1 0 - 6 1 0 - 7 1 0 - 8 1 0 - 9 U s n u m [mm⋅s-1] 1.7 1.8 1.9 2.0Figure 12
Slug velocity versus diffusivityD21 [m2⋅s−1] in logarithmic scale; solid line is represented by log⁡Us≈-0.023522log⁡D21+0.20714. The radius of the capillary tube is 0.4 [mm], and the successive lengths are l1=16 [mm] and l2=16.5 [mm].As a result, incrementing the diffusivity constantD21>0 by δD21>0 yields(38)γt,πD21+δD21<γt,πD21for  givenπ,t∈0,T.On the whole, augmenting D21>0 therefore compels the propulsion force, which is proportional to the surface tension difference γt,π, to decrease in magnitude. As a consequence, the transport coefficient D21>0 therefore does have bearing on the motion of the slug.In other words, augmenting the magnitude of the diffusivity constantD21>0 results in not only diminishing the magnitude of the slug velocity Us>0, but also causing the propulsion force to vanish/die earlier. As shown in Figure 12, one clearly sees that the π-dependent functional χπ now restricted to R,ls,μ=Consts evolves logarithmically with the diffusion constant D21>0 according to the following log-log relationship:(39)log⁡χπ≈-KD21log⁡D21+KD21,0whereKD21,KD21,0=Consts.As pointed out so far, it is understood that the scalar fields KD21, KD21,0 are partially π-dependant constants. This is an important feature which has not been considered in any of the two papers of Bico and Quéré [1, 5]. It may be remarked further that (39) tends to prove that the logarithm of the speed of the slug log⁡Us (since χπ≡Us) is proportional to the logarithm of the diffusivity constant log⁡D21. Moreover, since lt≈D21t (that diffusion length which provides a measure of how far the chemical concentration has propagated into the slug in the z-direction by diffusion in time t>0), the existence of a reference time τref>0 can be assumed such that(40)Ordlog⁡Us~Ordlog⁡dldt~log⁡D21τref+Const.,from which the shape of (39) at once comes to light. Physically speaking, if the diffusion constant D21>0 is to furnish us with a means of strengthening the self-propulsion mechanism, it is evident that it should below in magnitude. Hence, this establishes a relation between the diffusion constant D21>0 and the speed Us>0 of the slug. The next section considers the speed Us>0 of the slug versus the viscosity constant μ>0.


## 7.4. Case 2.4: Motion of Slug versus Viscosity

The nature of the movement of the slug with the character of the viscous force is here considered by establishing a relationship relating the speedUs>0 of the slug to the viscosity constant μ>0. To effect this purpose we ran several simulations, changing the magnitude of the property μ>0 at each run while restricting other property constants to those values appearing in Table 1. Further, at each run the pair of scalar fields μ,χπ is noted. In Figures 13 and 14 are plotted the pairs of scalar fields μi,χπi=1,2 on a log-log scale. As one can see from Figures 13 and 14, the dependence of the functional χπ finally restricted to R,ls,D21≡Consts on μi>0 (with i=1, 2) evolves on a log-log scale according to(41)log⁡χπ≈-Kμlog⁡μ+Kμ,0whereKμ,Kμ,0=Consts.As pointed out so far, it is understood that the scalar fields Kμ, Kμ,0 are partially π-dependant constants. Moreover, because this expression of Us>0 in terms of D21>0 is algebraically similar to that established in the foregoing section, our thoughts are thereby led to regard a decrease in μ>0 resulting in an exponential increase in Us>0, which manifests the nature of viscous effect.Figure 13
Slug velocity versus viscosityμ2 [mPa⋅s] in logarithmic scale; numerical result, approximated by log⁡Us≈-0.63452log⁡μ2+2.6509, is graphed in solid line; experimental result extracted from earlier work approximated by logUs≈-0.57214logμ2+2.39 is in dotted line. The radius of the capillary tube is 0.51 [mm], and the successive lengths are l1=16 [mm] and l2=16.5 [mm]; also, the viscosity μ1=17 [mPa⋅s] is kept constant.Figure 14
Slug velocity versus viscosityμ1 [mPa⋅s] in logarithmic scale; numerical result, approximated by logUs≈-0.42437log⁡μ1-2.0723, is graphed in solid line. The radius of the capillary tube is 0.51 [mm], and the successive lengths are l1=16 [mm] and l2=16.5 [mm]; also, the viscosity μ2=17 [mPa⋅s] is kept constant.To the informed reader the lawUs∝μ-1 is obvious whenever the Poiseuille law is pronounced in the mathematical development but, to the uninformed reader, the law Us∝μ-1/2 is unobvious. Thus, before advancing further, it might be interesting to justify this point.Since, in the earlier section (see Section7.2), we showed that Us∝ls-1 and, from Poiseuille flow, we have Us∝μ-1, it is plain that(42)Us∝1μ1α·μ21-αifμ=μ1,μ1=μ2;α∈0,1∝1μ1αifμ2=Const.;α∈0,1∝1μ21-αifμ1=Const.;α∈0,1.Equation (42) states that the dependence of the slug speed on the dynamic viscosity μi>0 of the i-fluid (with i=1, 2) is such that Us∝μi-β for some real β∈0,1. Moreover, it also states that if such dependence is known for one fluid, the other can be deduced, since whenever Us∝μ1-β necessarily Us∝μ21-β. Otherwise the condition Us∝μ-1 will be violated. Taking advantage of the properties of logarithms to convert the second and third lines of (42) into linear relationships, one obtains(43)log⁡Us=κi·log⁡μi+Const.∀i∈1,2,where∑i=1,2κi=1.Through numerical experiments, we deduced κ1,κ2=-0.42437,-0.63452. Since ∑i=1,2‍κi=1.05889, it follows that ∑i=1,2‍κi≈1, hence satisfying (43) and, therefore, (42). We have thus justified the results appearing in Figures 13 and 14.Having studied the motion of the slug under variations of the property vectorπ=(R,ls,D21,μ)T, we propose in the next section an algebraic expression, characteristics of the nature of the movement of the slug, on this ground.


## 8. Power-Law Slug Dynamics

In the foregoing sections, we examined the response of the slug subjected to the property constantπ. We also illustrated through case by case studies the effects of the latter on the movement of the slug. This section aims at proposing a power-law model, which might well be termed theLaw of Movement of the Slug. If we recall attention to the foregoing sections as regards the nature of the functional χπ with respect to the property π, we shall see that the speed Us>0 of the slug is related to the elements of the property vector π=(R,ls,D21,μ)T as(44)Us≈KRR+KR,0,KR,KR,0=Consts;see35;-Klslog⁡ls+Kls,0,Kls,Kls,0=Consts;see36;-KD21log⁡D21+KD21,0KD21,KD21,0=Consts;see39;-Kμlog⁡μ+Kμ,0Kμ,Kμ,0=Consts;see41.On the other hand, from the above graphical results it is observed that the quantities KR, Kls, |K, and Kμ, respectively, roughly interrelate themselves as follows:(45)OrdKR~OrdKls~Ord2Kμ~Ord50KD21~1.If we pause ourselves for a moment and ponder on the above relationship, we would obviously conclude that the relative effect arising from the diffusion constant D21>0 is little. Thence, paralleling (31) with that which one would deduce at leading order in a fully developed slug flow in a circular pipe of the same radius R together with the neglect of log⁡D21 yields(46)Us~ξRγ1,0-γ2,04μlsα,Ueff~ξRγ1,0-γ2,04μlsα-1,where α=Kls~Kμ>0 and ξ>0. Writing the speed Us>0 of the slug as such puts forward the term of the concept of aPoiseuille flow. In these, the scalar fields ξ>0, α define theslug dynamics consistency andslug behaviour indices (with α≈1.0); the scalar field Ueff>0 is aneffective slug velocity. To deduce the scalar field ξ>0 on this basis one may then proceed through practical experimentation.In the paper of Lunati and Or [35], the authors mentioned that the slug velocity, for the case of a gravity-driven flow in a capillary tube, is lower (without any precise order of magnitude) than deduced in a Poiseuille flow. Here, the result points to a contrary conclusion, for using the data given in Table 1 one demonstrates that the speed Us>0 of the slug isgreater by an ϵ-order ϵ~10-2, roughly. Such discrepancy might be due to the neglect of dissipative effects occurring at the wedges of the slug.To this end, it seems therefore plausible to assume the following general description:(47)Xst,π≈Xs,∞1-exp⁡-ξt-t0Xs,∞Rγ1,0-γ2,04μlsα.We agree here that the property π includes all property constants appearing on the right-hand of (47), as can be easily understood. There is an immediate consequence which follows from (47). Setting λ=Rγ1,0-γ2,0/4μls, one concludes that ∂tXs∝λ in the early-time (linear) regime, whereas ∂tXs∝exp-λ in the late-time (nonlinear) regime. Since a Poiseuille flow is one such that its flow fields satisfy u∝λ, it follows that the slug flow obeys Poiseuille flow only in the early-time regime.To round off the present discussion, we point out that, withξls=30[mm]≈22.3633 and ξls=60[mm]≈13.2664 along with the numerical data given in Tables 2 and 3, plotting of (47) coincides with those given (in solid lines) in Figure 15, for comparing (29) with (47) gives(48)ξπ≈η1η2Rγ1,0-γ2,04μls1-α≈1for  givenπ,sinceα≈1.Finally, combining (23) with (47) and (48), one arrives at(49)∂tXs≈Us1-XsXs,∞≤Us,dXsdγ≈Ueff1-XsXs,∞≤Ueff.Hence, it would appear that (49) hold for every fixed instant of time t>0, provided that ∂tγ≈Const. Hence, if we let the scalar field Up>0 be the Poiseuille velocity field, one deduces(50)∂tXs≈Rγ1,0-γ2,04μlsdXsdγ≈UpdXsdγwhen∂tγ≈Const.,since Up≈Us·Ueff-1. As pointed out before, it then follows that OrddXs/dγ~ϵ-2.Figure 15
Results with (47) and numerical data. Equation (47) is described in solid line for ls=30 [mm] and in dotted line for ls=60 [mm]. Likewise, numerical data is described with circles for ls=30 [mm] and with squares for ls=60 [mm].In Figure16 are displayed the flow contours generated by the movement of the slug into its capillary tube for t=5k/2, where k∈1,4, the time evolution of the advancing interface in the early-time regime. The region Ωs=Ω1∪Ω2 designates the domain partially filled with the miscible fluid; the region Ωv+ designates that occupied by the vapour phase fluid.Figure 16
Illustration of the flow contours at the sequence of timestkk∈1,5=5k/2k∈0,4 [s].In Figure17 are displayed the distribution of the fraction volume of fluids ϕ12r,t within the regions Ω1∪Ω2 and Ωv-∪Ωv+, the time evolution of the advancing interface in the early-time regime, for the sequence of time tk=2k, where k∈0,5. This completes the discussion of this part.Figure 17
Evolution of the fraction volume of fluid to illustrate the front interfaceΓ1 of the slug at the sequence of times tkk∈1,5=2kk∈0,5 [s].


## 9. Conclusion

Hitherto, it is believed insofar as we are knowledgeable that this is the first time that such a study has been undertaken. Moreover, perusing our results, one admittedly concludes that ours are very similar to the observations reported by Bico and Quéré [1, 5].On the whole, we have observed the following points.(i)
A mathematical model is developed and implemented in COMSOL Multiphysics.(ii)
Two regimes were identified: the early-time (linear) and late-time (nonlinear) regimes. In the early-time regime, the speed of the slug is proportional to the Poiseuille flow, whereas in the late-time regime it is proportional to the exponential of the Poiseuille flow.(iii)
Exploring the linear regime parametrically, the dependence of the speed of the slug on its sizels, diffusivity D21 and dynamic viscosity μ constants, and the radius R of the capillary tube is found to be power-law functions. Moreover, the speed is found to be insensitive to D21. These parametric studies show that the size and properties of the slug do have bearing on its motion.(iv)
For an initially chosen set of physical properties, the position of the slug evolves exponentially with respect to time according to (47).Of the above statements, the second one confirms the conclusion reached by Piroird et al. [43] in one of their papers, namely, that on a long time scale two regimes were observed; and the third one confirms that found in another paper of Piroird et al. [44], namely, that by varying the geometrical properties associated with the droplet and its capillary tube, the nature of the flow changes.Additionally, it should be stressed that it is notnecessary to argue or even to enumerate in the context of practical applications, for we believe that the above mathematical development issufficient to prove its usefulness in practical applications. However, it should be borne in mind that, without the numerical data ξ, α, and Xs,∞, the proposed model is insufficient to allow one to predict the position of the slug at some time t0≥0 thereafter. Nonetheless, it does indicate that in this range it is not possible to describe the dynamics of the slug without knowing, for instance,at what time,t∞>0say, the slug must reach its destination,Xs,∞>0say. Such is the conclusion to which we are led.Taking into account the foregoing discussions, it may be declared that our model, namely, (47), can do good service in the selection of microfluidic (or more preciselycapillary) devices. To conclude, it may be said that the outcome of the present part results in the establishment of the following proposition.Proposition 1.
Givenκ1,κ2, there is a pair of π-dependent property constants. By a proper choice of the property π=π1,…,πd∈Rd, the path of a self-propelling slug in a capillary tube obeys the following law:(51)Xπ;t=κ11-exp⁡-κ2twithκi=fiπ,fori=12.We name (51) thePrototype Slug Self-Propulsion Model Equation because it gave agreement to both time and physical properties of the slug.



---

*Source: 1234642-2016-10-27.xml*