# Scaling Laws of Droplet Coalescence: Theory and Numerical Simulation

Authors: M. Irshad Khodabocus; Mathieu Sellier; Volker Nock
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

When two Newtonian liquid droplets are brought into contact on a solid substrate, a highly curved meniscus neck is established between the two which transforms the bihemispherically shaped fluid domain to a hemispherically shaped domain. The rate at which such topological transformation, called coalescence phenomenon, evolves results from a competition between the inertial force which resists the transformation, the interfacial force which promotes the rate, and the viscous force which arrests it. Depending on the behaviour of these forces, different scaling laws describing the neck growth can be observed, predicted theoretically, and proved numerically. The twofold objective of the present contribution is to propose a simple theoretical framework which leads to an Ordinary Differential Equation, the solution of which predicts the different scaling laws in various limits, and to validate these theoretical predictions numerically by modelling the phenomenon in the commercial Finite Element software COMSOL Multiphysics.


---

## Body

## 1. Introduction

When two Newtonian liquid droplets are brought into contact on a solid substrate, a highly curved meniscus bridge is established and thereafter driven by surface tension so as to transform the bulk flow region from bihemispherically shaped at early-time regimes to more hemispherically shaped at late-time regimes. Such topological transformation, involving the merging of two or more sessile droplets into a single one on a solid surface, is calleddroplet coalescence phenomenon.Physically speaking, the coalescence process consists of three regimes, called theviscous,visco-inertia, andinertial coalescence regimes [1]. In each regime, a range of surprisingly complex behaviours can be observed, so to speak topological changes of the interface [2]. Thus, the subject of droplets coalescence is a fascinating, multifaceted subject of inquiry with many interesting phenomena to learn and discover. Indeed, the subject has intrigued many physicists [3–6] and mathematicians [7, 8].In actual fact, authors of [9] appear to be the first scientists who have addressed the subject. The authors study the formation of vortex rings by droplets falling into liquids and some allied phenomena. For a miscible droplet falling from not too great a height into its bulk fluid, it is observed that the droplet descends through the bulk fluid as a ring, whereas for the case of an immiscible droplet, [9] observed the droplet to descend through the bulk fluid as a spherical droplet. Thereafter, there have been many attempts on experimental, theoretical, and numerical sides to elucidate the underlying physics governing the dynamics of coalescence.Authors of [10] have studied the partial coalescence of a water droplet, with radius R∈0.02,1.3 mm, with the free-surface of a bulk fluid composed of water-glycerol mixtures experimentally. For cases R→0.02, 1.3 mm, the authors drew attention that the coalescence process is arrested by gravitational and viscous effects and, depending on these two effects, three coalescence regimes were observed, which [10] termedgravity,inertio-capillary, andviscous regimes. Furthermore, the authors argued that the coalescence time is strongly dependent upon the Bond number Bo=4ρ1-ρ2gR2/γ and the Ohnesorge number Oh=μ1/2ρ1+ρ2γR1/2, where ρ1, μ1, and γ denote the density, dynamic viscosity, and surface tension of the droplet, respectively; ρ1 and μ1 denote the properties of its vapour phase surrounding, and g denotes the gravitational acceleration.Reference [7] addressed the more general case of the coalescence of droplets surrounded by a viscous outer fluid, both analytically and numerically, using asymptotic methods and the integral equation given by [11]. Describing the initial coalescence regime by the Stokes equations, they showed that the three-dimensional solution is asymptotically equivalent to the two-dimensional one reported in the works of [12–16]. When the viscosity of the outer fluid is disregarded, the radius and width of the highly curved meniscus bridge were found to obey the following laws: rm~tγ/πμln⁡[tγ/μR] and Δ∝rm3, respectively, where t is the time, and properties R, γ, and μ are the radius, surface tension, and viscosity constant of the droplet; when the viscosity of the outer fluid is taken into account, a toroidal bubble of radius Δ∝rm3/2 is formed by the outer fluid, and rm~tγ/4πμln⁡[tγ/μR]. Employing the Euler equations in the inertial regime, they showed that rm∝γR/ρ1/4t1/2, where ρ denotes the density of the coalescing droplets. (In the present study, the property bπ;t will stand for rm; π designates the set of property constants upon which it depends.)The paper of [7] is one classic paper which treats the coalescence of droplets both theoretically and numerically. The method used by [7] is termedconventional by [8]. The work of [8] is another classic paper on the subject. They addressed the subject using the conventional method and a different method calledinterface formation/disappearance method and comparing their corresponding results with experimental findings. The authors argued that their proposed model shows better agreement with experimental data than the conventional one.In the literature, [4] has given a combined theoretical and experimental account of the coalescence of water droplets on a solid substrate surrounded by an atmosphere of nitrogen saturated with water. For small contact angle θ<20∘ and large contact angle θ∈70∘,90∘, the bidroplet undergoes very fast coalescence. Moreover, the relaxation time characterising the coalescence flow is found to be tc∝R (or, according to [4], tc=πμ/3Kγ), and the bulk capillary relaxation time, tb=Rμ/γ≈tc×10-7, where R, γ, and μ are the radius, surface tension, and viscosity constant of the bidroplet system.From a hydrodynamical viewpoint, [3] found that the time evolution of the characteristic scale bπ;t of that small liquid bridge—where the meniscus of the bidroplet system is highly curved—obeys a time-dependent scaling/power law: bπ;t∝tk, where k>0. For instance, in theinertial coalescence regime it was found that k≈1/2, whereas in theviscous coalescence regimek≈1. The confirmation of these laws as regards the growth of the property bπ;t in those coalescence regimes was also reported by earlier scientists. Reference [2] addressed the early-time evolution of the highly curved meniscus bridge experimentally. They found that bπ;t∝K1/4t1/2, where K is the inverse of the curvature of the meniscus bridge; by comparison, K=γR/ρ [7].On increasing the viscosity of the fluid, [17] has shown the linear dependence of bπ;t on time. Reference [18] has studied the influence of the geometrical properties of low viscosity liquid droplets experimentally. For contact angles θ<90∘, they deduced bπ;t∝t2/3, while for θ>90∘ their analysis yielded to bπ;t∝t1/2. Most interestingly, they demonstrated the 1/2 exponent scaling law in the limit that θ→90∘, thereby unifying the coalescence of liquid droplets and freely suspended droplets in the inertial regime. The 2/3 exponent scaling law was also argued by [19] who attributed it to the self-similar description of the shape of inviscid pinchoff of a liquid droplet in a still air.On the other hand, [20] has considered that the influence of surface tension gradient on droplet coalescence was addressed on a combined experimental and numerical basis. The situation considered by these authors was the coalescence of a droplet deposited on the surface of a miscible liquid reservoir. Their results reveal three distinct coalescence regimes based on the reservoir-to-droplet surface tension ratio λγ>0. For instance, the situation λγ<0.42 results in the ejection of a small daughter droplet from the top of the top of the coalescing droplet. Contrarily, it was only when λγ∈0.42,0.93 that total coalescence could be observed.Reference [21] has addressed the early-time coalescence of viscous droplets on a flat, wettable substrate. The authors have shown in terms of their respective radii R0, initial heights h0≪R0, surface tensions γ>0, and dynamics viscosities μ>0 that bπ;t≈γh03/μR021/2t, which led them to conclude that the evolution of bπ;t is highly sensible to the geometrical properties of the liquid droplets. Reference [22] considered the coalescence of air bubbles and salt water droplets systems in silicone oils experimentally. The proportionality constant and the α-exponent of the scaling law bπ;t∝tα are analysed; the authors argued that the evolution of bπ;t is independent of the viscosity of its outer fluid, changing only the proportionality constant by a factor roughly equal to 1.5 and leaving the α-exponent unaffected.The above references reveal extensive studies of the characters of the propertybπ;t=κ·tα with respect to some properly/carefully chosen pair of scalar fields κ,α, paying little attention to the characters of the coalescence flow fields in three-dimensional space and on cutting planes and leaving the simultaneous establishment of these power law regimes from one and only one generic equation untouched. Moreover, should an attempt be made to study a class of physical situations ranging from droplets spreading to self-propulsion and coalescence, the several modelling frameworks proposed by those scientists are inconvenient to employ without reconstructing either the mathematical framework or the numerical framework, though the outcomes of their works are no doubt striking. Thus, on the one hand, we must admit that a b-evolution equation which permits the deduction of those power laws has never been derived and, on the other hand, the models proposed by those authors make it difficult, if not impossible, to display the global characters of the coalescence flow fields as well as to adapt from one physical situation to another.In this part, we develop a sound framework which upon variations of an interdependent leading parameter permits to study either the coalescence of two liquid droplets or the spreading of a single droplet on a three-dimensional solid substrate. Taking advantage of the proposed geometrical model, we then establish through theoretical analysis thatb-evolution equation and derive from it the power law relationships bπ;t∝tα (α=1/2, 2/3, 1).To demonstrate that the proposed model might be adopted with advantage, the1/2, 2/3, 1-exponent scaling laws are proved numerically, and the characters of the coalescence flow fields are illustrated in three-dimensional space and on cutting planes.


## 2. Theory

### 2.1. Mathematical Formulations

To formulate the problem mathematically, we reason in the following manner. Suppose two liquid dropletsΩ1, Ω2⊂R3 coalesce on a solid substrate S⊂R3, at some time t>0 they will occupy amaster domain ⋃i=1,2‍Ωi satisfying the geometrical condition(1)vol⋃i=1,2Ωi≤volΩ1+volΩ2forallt∈0,T.Following [8], equality follows from (1) only when the size of theinterfacial circular sector centered at the point xs∈S (see Figure 1) is zero, and strict inequality whenever it is nonzero, as is easily understood. In the present work, the domain ⋃i=1,2‍Ωi, which can vary in time, will be termedcoalescence domain, though the termbidroplet system is also attributed to; the set 0,T (with 0<T<∞) denotes the period of time over which the study is carried out. Mathematically, the coalescence domain is bounded by two disjoint sets called theinterface∂(⋃i=1,2‍Ωi) and thefootprint⋃i=1,2‍Ωi∩S of the coalescence domain. In Figure 1, we illustrate the coalescence of two viscous liquid droplets Ω1, Ω2⊂R3 conditioned byzero-flux of momentum, mass and energy across the (free) surface ∂(⋃i=1,2‍Ωi) and the (contact) surface ⋃i=1,2‍Ωi∩S; the null flux vectors φi=0 (with i=1, 2, s) stand for those boundary conditions. Figure 2 exposes to view the projection of the coalescence domain in the y,z-plane (top) and in the x,z-plane (bottom). Points x0, xs∈S designate the rear and centroid of the bidroplet system; and, points x1, x2∈S the centroids of its first and second droplets, both taken at the instant of time t=0. Quantities R1, R2 (with R2≤R1) are the initial radii of the coalescing droplets, and b0≡bπ;t=0 the characteristic size at their corresponding points of contact. In particular, it is assumed in Figure 2 that the coalescence domain is initiallybihemispherical, L units in length in the direction of ey, with respective radii R1, R2 (with R2≤R1, by hypothesis), and satisfying the condition R1+R2<L/2. (We agree that R1+R2=L/2 necessarily implies that b0=0 and, R1+R2>L/2 implies that the coalescing droplets have not yet been brought into contact; see Figure 2.)Figure 1
A freehanded illustration of the coalescence of two liquid droplets.Figure 2
The projection of the coalescence domain in they,z-plane (top) and in the x,z-plane (bottom).The dynamic state of the coalescence domain is characterised completely by the following macroscopic vector and scalar fields quantities: thevelocityux,t=ux,uy,uzT and thepressurepx,t. We draw attention that, the vector field ux,t describes the coalescence flow within the coalescence domain ⋃i=1,2‍Ωi.Thegoverning equations describing these property variables in ⋃i=1,2‍Ωi are the following generic momentum and mass transport equations of a divergence-free flow field:(2)∂tρu+divρu⊗u+∇p=ρfγ+divSuin⋃i=1,2Ωi×0,T;divρu=0in⋃i=1,2Ωi×0,T.In the above system of balance laws, the scalar field px,t denotes the pressures field and the vector field fγx,t denotes the density of volume force due to the surface tension. The tensor Su denotes theincompressible Newtonian viscous stress tensor, defined as(3)Su=μgradu+graduTin⋃i=1,2Ωi×0,T.In (2), the property constants ρ, μ>0 are, respectively, thedensity anddynamic viscosity of the mixed fluid occupying the domain ⋃i=1,2‍Ωi. The operator ⊗ designates thetensor product operator; for given vector field variables u, v∈R3, their tensor product yields to a 3×3 matrix with the following entries: u⊗v=uivj (with i, j=x,y,z).To model the vector fieldfγx,t, we employ the interfacial scalar field γx,t into the following alternative expression(4)fγx,t=divγ∇s·nn-∇sγon∂⋃i=1,2Ωi×0,T;∇s·=I-n⊗n∇·.The detailed derivation of (4) is found in [23]. The vector nx,t=nx,ny,nzT is the unit outer vector normal to ⋃i=1,2‍Ωi and ⋃i=1,2‍Ωi∩S (≠⋃i=1,2‍Ωi, since, Ωi∩S≠Ωi for i=1, 2), respectively. The operator ∇s· is thesurface gradient operator on ∂(⋃i=1,2‍Ωi), and the matrix I stands for theidentity matrix.


### 2.2. Boundary and Initial Conditions

At the instant of timet=0 the distributions of the velocity field ux,t, the pressure field px,t and the chemical concentration cx,t should be supplemented in ⋃i=1,2‍Ωi, and thereafter conditioned on boundaries ∂(⋃i=1,2‍Ωi) and ⋃i=1,2‍Ωi∩S for all t∈0,T. These are discussed in the following two sections.The dynamic boundary conditions imposed on the interface∂(⋃i=1,2‍Ωi) for the velocity field ux,t are the Neumann-type conditions(5)Su-pIn=-pextI+γ∇s·nInon∂⋃i=1,2Ωi×0,T,Su-pIt=∇sγon∂⋃i=1,2Ωi×0,T.In the above conditions, the property pext>0 designates theexternal pressure applied on the interface ∂Ω1∪Ω2; γ∇s·n is the force per unit area due to local curvature of that interface. The term on the right-hand side of the second equation designates the tangential stress associated with gradients in the surface tension coefficient. The proper boundary conditions to be considered on the contact surface ⋃i=1,2‍Ωi∩S are now discussed. At some time t>0 thereafter, the footprint of the coalescence domain ⋃i=1,2‍Ωi is expected to deform from abidisc shaped domain to amonodisc shaped domain. Therefore, a wetted wall boundary condition is considered, thesolid-fluid boundary conditions. Thus, if the scalar quantity β>0 designates theNavier-slip coefficient [24], considering theNavier-slip-with-friction boundary condition yields(6)Sun+βu·t=0on⋃i=1,2Ωi∩S×0,T,u·n=0on⋃i=1,2Ωi∩S×0,T.In (5), (6), the vector tx,t=tx,ty,tzT designates the tangent vectors to boundaries ∂(⋃i=1,2‍Ωi) and ⋃i=1,2‍Ωi∩S. The following initial condition is imposed on the vector field ux,t at t=0:(7)ux,t=u0xin⋃i=1,2Ωi×0,The momentum and mass transport equations along with the foregoing set of boundary and initial conditions are solved using the COMSOL Moving Mesh Interface, discussed in what follows.


### 2.3. Moving Coalescence Domain Explained

At this stage, we have beforehand the mathematical framework and its respective domain of definition. Bearing in mind that the bulk motion of the coalescence domain results from topological changes of its interface from one instant of time to another, to adapt the mathematical descriptions with the movement of the coalescence domain⋃i=1,2‍Ωi we employ themoving mesh method, which we shall now throw some light upon. For more details, see the works of [23, 25, 26].Let the interface∂(⋃i=1,2‍Ωi) of the coalescence domain designates amaterial point-set consisting of material points X=X,Y,ZT at the instant of time t=0. Under this condition, let this coalescence domain at some time t≥0 be described as aspatial point-set consisting of points x=x,y,zT. On the other hand, the coalescence domain being subjected to smooth deformation from one instant to another, there follows that(8)x=xX,t,on∂⋃i=1,2Ωit=0×0,T,where X=xX,0. Having thesematerial andspatial vector fields at hand, to relate the topological changes of the interface to its bounded bulk fluid, it suffices to equate the associated velocities on the said interface. Thus, if the vector field Ux,t=Ux,Uy,UzT describes the movement of the interface, there results in(9)u·n=U·n,on∂⋃i=1,2Ωi×0,T;x=xX,t.Note that, we write ∂(⋃i=1,2‍Ωit=0) instead of ∂(⋃i=1,2‍Ωi) to arrest attention that the interface is undeformed initially.


## 2.1. Mathematical Formulations

To formulate the problem mathematically, we reason in the following manner. Suppose two liquid dropletsΩ1, Ω2⊂R3 coalesce on a solid substrate S⊂R3, at some time t>0 they will occupy amaster domain ⋃i=1,2‍Ωi satisfying the geometrical condition(1)vol⋃i=1,2Ωi≤volΩ1+volΩ2forallt∈0,T.Following [8], equality follows from (1) only when the size of theinterfacial circular sector centered at the point xs∈S (see Figure 1) is zero, and strict inequality whenever it is nonzero, as is easily understood. In the present work, the domain ⋃i=1,2‍Ωi, which can vary in time, will be termedcoalescence domain, though the termbidroplet system is also attributed to; the set 0,T (with 0<T<∞) denotes the period of time over which the study is carried out. Mathematically, the coalescence domain is bounded by two disjoint sets called theinterface∂(⋃i=1,2‍Ωi) and thefootprint⋃i=1,2‍Ωi∩S of the coalescence domain. In Figure 1, we illustrate the coalescence of two viscous liquid droplets Ω1, Ω2⊂R3 conditioned byzero-flux of momentum, mass and energy across the (free) surface ∂(⋃i=1,2‍Ωi) and the (contact) surface ⋃i=1,2‍Ωi∩S; the null flux vectors φi=0 (with i=1, 2, s) stand for those boundary conditions. Figure 2 exposes to view the projection of the coalescence domain in the y,z-plane (top) and in the x,z-plane (bottom). Points x0, xs∈S designate the rear and centroid of the bidroplet system; and, points x1, x2∈S the centroids of its first and second droplets, both taken at the instant of time t=0. Quantities R1, R2 (with R2≤R1) are the initial radii of the coalescing droplets, and b0≡bπ;t=0 the characteristic size at their corresponding points of contact. In particular, it is assumed in Figure 2 that the coalescence domain is initiallybihemispherical, L units in length in the direction of ey, with respective radii R1, R2 (with R2≤R1, by hypothesis), and satisfying the condition R1+R2<L/2. (We agree that R1+R2=L/2 necessarily implies that b0=0 and, R1+R2>L/2 implies that the coalescing droplets have not yet been brought into contact; see Figure 2.)Figure 1
A freehanded illustration of the coalescence of two liquid droplets.Figure 2
The projection of the coalescence domain in they,z-plane (top) and in the x,z-plane (bottom).The dynamic state of the coalescence domain is characterised completely by the following macroscopic vector and scalar fields quantities: thevelocityux,t=ux,uy,uzT and thepressurepx,t. We draw attention that, the vector field ux,t describes the coalescence flow within the coalescence domain ⋃i=1,2‍Ωi.Thegoverning equations describing these property variables in ⋃i=1,2‍Ωi are the following generic momentum and mass transport equations of a divergence-free flow field:(2)∂tρu+divρu⊗u+∇p=ρfγ+divSuin⋃i=1,2Ωi×0,T;divρu=0in⋃i=1,2Ωi×0,T.In the above system of balance laws, the scalar field px,t denotes the pressures field and the vector field fγx,t denotes the density of volume force due to the surface tension. The tensor Su denotes theincompressible Newtonian viscous stress tensor, defined as(3)Su=μgradu+graduTin⋃i=1,2Ωi×0,T.In (2), the property constants ρ, μ>0 are, respectively, thedensity anddynamic viscosity of the mixed fluid occupying the domain ⋃i=1,2‍Ωi. The operator ⊗ designates thetensor product operator; for given vector field variables u, v∈R3, their tensor product yields to a 3×3 matrix with the following entries: u⊗v=uivj (with i, j=x,y,z).To model the vector fieldfγx,t, we employ the interfacial scalar field γx,t into the following alternative expression(4)fγx,t=divγ∇s·nn-∇sγon∂⋃i=1,2Ωi×0,T;∇s·=I-n⊗n∇·.The detailed derivation of (4) is found in [23]. The vector nx,t=nx,ny,nzT is the unit outer vector normal to ⋃i=1,2‍Ωi and ⋃i=1,2‍Ωi∩S (≠⋃i=1,2‍Ωi, since, Ωi∩S≠Ωi for i=1, 2), respectively. The operator ∇s· is thesurface gradient operator on ∂(⋃i=1,2‍Ωi), and the matrix I stands for theidentity matrix.


## 2.2. Boundary and Initial Conditions

At the instant of timet=0 the distributions of the velocity field ux,t, the pressure field px,t and the chemical concentration cx,t should be supplemented in ⋃i=1,2‍Ωi, and thereafter conditioned on boundaries ∂(⋃i=1,2‍Ωi) and ⋃i=1,2‍Ωi∩S for all t∈0,T. These are discussed in the following two sections.The dynamic boundary conditions imposed on the interface∂(⋃i=1,2‍Ωi) for the velocity field ux,t are the Neumann-type conditions(5)Su-pIn=-pextI+γ∇s·nInon∂⋃i=1,2Ωi×0,T,Su-pIt=∇sγon∂⋃i=1,2Ωi×0,T.In the above conditions, the property pext>0 designates theexternal pressure applied on the interface ∂Ω1∪Ω2; γ∇s·n is the force per unit area due to local curvature of that interface. The term on the right-hand side of the second equation designates the tangential stress associated with gradients in the surface tension coefficient. The proper boundary conditions to be considered on the contact surface ⋃i=1,2‍Ωi∩S are now discussed. At some time t>0 thereafter, the footprint of the coalescence domain ⋃i=1,2‍Ωi is expected to deform from abidisc shaped domain to amonodisc shaped domain. Therefore, a wetted wall boundary condition is considered, thesolid-fluid boundary conditions. Thus, if the scalar quantity β>0 designates theNavier-slip coefficient [24], considering theNavier-slip-with-friction boundary condition yields(6)Sun+βu·t=0on⋃i=1,2Ωi∩S×0,T,u·n=0on⋃i=1,2Ωi∩S×0,T.In (5), (6), the vector tx,t=tx,ty,tzT designates the tangent vectors to boundaries ∂(⋃i=1,2‍Ωi) and ⋃i=1,2‍Ωi∩S. The following initial condition is imposed on the vector field ux,t at t=0:(7)ux,t=u0xin⋃i=1,2Ωi×0,The momentum and mass transport equations along with the foregoing set of boundary and initial conditions are solved using the COMSOL Moving Mesh Interface, discussed in what follows.


## 2.3. Moving Coalescence Domain Explained

At this stage, we have beforehand the mathematical framework and its respective domain of definition. Bearing in mind that the bulk motion of the coalescence domain results from topological changes of its interface from one instant of time to another, to adapt the mathematical descriptions with the movement of the coalescence domain⋃i=1,2‍Ωi we employ themoving mesh method, which we shall now throw some light upon. For more details, see the works of [23, 25, 26].Let the interface∂(⋃i=1,2‍Ωi) of the coalescence domain designates amaterial point-set consisting of material points X=X,Y,ZT at the instant of time t=0. Under this condition, let this coalescence domain at some time t≥0 be described as aspatial point-set consisting of points x=x,y,zT. On the other hand, the coalescence domain being subjected to smooth deformation from one instant to another, there follows that(8)x=xX,t,on∂⋃i=1,2Ωit=0×0,T,where X=xX,0. Having thesematerial andspatial vector fields at hand, to relate the topological changes of the interface to its bounded bulk fluid, it suffices to equate the associated velocities on the said interface. Thus, if the vector field Ux,t=Ux,Uy,UzT describes the movement of the interface, there results in(9)u·n=U·n,on∂⋃i=1,2Ωi×0,T;x=xX,t.Note that, we write ∂(⋃i=1,2‍Ωit=0) instead of ∂(⋃i=1,2‍Ωi) to arrest attention that the interface is undeformed initially.


## 3. Numerical Implementation

In this section, the construction of the physical domain⋃i=1,2‍Ωi, described in the x,y,z-space by Figure 1 and in the y,z-plane by Figure 2, and the numerical implementation of those equations and conditions, namely, (1) to (9), describing its dynamic state is outlined.The bidroplet system—that is, the coalescence domain—is built as a three-dimensional bihemispherical computational domain sitting on ax,y-plane solid substrate. The solid substrate is defined as a rectangular x,y-plane sized Li>0 units in length in the direction of ei (with i=x, y). To realise these geometrical constructions, the radii R1, R2, and position vectors xk=xk,yk,zkT∈⋃j=1,2‍Ωj∩S (with k=1, 2, s) must be supplied, and are tabulated in Table 1. They are based on the following line of reasoning.Table 1
Initial and boundary conditions imposed on the coalescence domains⋃i=1,2‍Ωi. Geometrical properties used to generate the b0-dependent coalescence domains ⋃i=1,2‍Ωi.
ρ μ γ p e x t u 0 x  [g⋅cm−3] [Pa⋅s] [N⋅m−1] [Pa] [m⋅s−1] 1.3 1.4 6 × 1 0 - 2 1.0 0,0 , 0 T Lx=2Ly R 1 = R 2 x 1 x 2 x s  [mm] [mm] [mm] [mm] [mm] 10 1.0 6,5 , 0 T 4,5,0T 5,5,0TInitially, the bidroplet system is assumed symmetrical with respect to they,z-plane. On the solid substrate, the position vector xs=0,ys,0T denotes its centroid of the coalescence domain ⋃i=1,2‍Ωi; the position vectors x1=0,y1,0T, x2=0,y2,0T denote the centroids of its first and second droplets, respectively. These descriptions are clearly represented in Figure 2. To stress the dependence of the property bπ;t on those geometrical quantities, we consider the following expressions:(10)bπ;t=b0=Ri2-yi-ys21/2att=0,where i=1, 2. Furthermore, to stress that the property b0>0 is a relatively small quantity initially, we impose the condition b0≪minR1,R2; in regard to the freehand sketches of Figure 2, this means that b0≪R2≤R1. We note that bπ;t≥b0 for all t≥0; thus the coalescence phenomenon is arelative process, since it is studied through the evolution of the relative property bπ;t-b0. Consequently, the expansion rate at that region where the curvature is very large is very dependent upon the magnitude of b0>0. Clearly, the smaller be the size of the property b0, the higher be the curvature at the common point of contact of two droplets, and, consequently, the faster be the coalescence rate at early-time.Topologically, the case for whichb0=0 prescribes aone-point interface, while for the other cases for which b0≠0, all prescribeinfinitely many-points interface. Since those interfaces disappear for all t∈0,T, the first case might well be termed theone-point model and the other three casesinfinitely many-points model.Granted the above descriptions, the following set of equations is implemented: the COMSOLLaminar Two-Phase Flow Moving Mesh equations set, having appropriate interfaces of their initial and boundary conditions. This model, as its name implies, is used to describe the spatiotemporal descriptions of the variable field function ux,t. To adapt the COMSOL u-equation model with the physical situation under consideration, the COMSOL default model in question is expanded so as to suit the u-set of initial and boundary conditions given in earlier (theoretical) sections. The values attributed to those conditions are reported in Table 1.A series of CFD experiments are ran on this ground and their results are presented in what follows.


## 4. Results and Discussion

The results to which the proposed model led on a combined theoretical and numerical basis are discussed here.


### 4.1. Preliminaries

The highly curved meniscus is a fundamental entity. On the grounds of earlier works [3, 18, 19, 21, 27], it can be characterised as follows. By hypothesis, let the d-dimensional sequence π=π1,…,πd denotes the property-set of d≥1 positive property constants, upon which rests the character of the coalescence dynamics. Further, let ψαπ be a π-dependent property, where the subscript α stands for the associated α-exponent law. Then, following earlier works one can postulate the following expression:(11)bπ,t=ψαπ·tα,fort∈0,T,whereα>0.Given L, T the fundamental units of length and time, the dimensional unit of the property ψαπ is deduced as follows: We observe that [∂tb]=[L·T-1], therefore, [ψα]=[L·T-α], where α>0. Moreover, if by properties b0, Ψref, τref>0 are meant the characteristic scaling of the scalar fields bπ;t, ψπ, t, the dimensionless form of (11) reads: bπ;t=tα, where Ordb0~Ψref·τrefα for some α>0, meaning that the law is independent of the working Newtonian fluids. Consequently, sufficient is that our CFD framework is capable of proving coalescence laws. These statements motivate us to study the coalescence of the bihemispherically shaped droplet system.


## 4.1. Preliminaries

The highly curved meniscus is a fundamental entity. On the grounds of earlier works [3, 18, 19, 21, 27], it can be characterised as follows. By hypothesis, let the d-dimensional sequence π=π1,…,πd denotes the property-set of d≥1 positive property constants, upon which rests the character of the coalescence dynamics. Further, let ψαπ be a π-dependent property, where the subscript α stands for the associated α-exponent law. Then, following earlier works one can postulate the following expression:(11)bπ,t=ψαπ·tα,fort∈0,T,whereα>0.Given L, T the fundamental units of length and time, the dimensional unit of the property ψαπ is deduced as follows: We observe that [∂tb]=[L·T-1], therefore, [ψα]=[L·T-α], where α>0. Moreover, if by properties b0, Ψref, τref>0 are meant the characteristic scaling of the scalar fields bπ;t, ψπ, t, the dimensionless form of (11) reads: bπ;t=tα, where Ordb0~Ψref·τrefα for some α>0, meaning that the law is independent of the working Newtonian fluids. Consequently, sufficient is that our CFD framework is capable of proving coalescence laws. These statements motivate us to study the coalescence of the bihemispherically shaped droplet system.


## 5. The Laws ofbπ;t Predicted Theoretically

Granted the proposed three-dimensional coalescence domain (see Figure1), and the accompanying two-dimensional geometrical model (see Figure 2) wherein are reflected its geometrical properties in the x,z-plane and in the y,z-plane, one can predict the time evolution of the property bπ;t in the power law regimes. In the present section, the theoretical prediction of the property bπ;t in those regimes is undertaken. Let it be given ϵ≪1, b0>0, and xs=xs,ys,zsT∈⋃i=1,2‍Ωi∩S, and consider the elementary domain defined by(12)Bϵ0xs=x,y,z:y-ys≤ϵ,x-xs+z-zs≤b0,with boundary ∂Bϵ0xs, at the instant of time t=0 (see Figure 3 for geometrical insight). Note that, the elementary domain Bϵ0xs⊂⋃i=1,2‍Ωi may be conceived of as a 180 degree sector of ϵ-order thick aligned with the x,z-plane, and having its footprint centered at the point xs. The region Bϵ0xs is characterised geometrically by the following length and width at time t=0: leng∂Bϵ0=πb0, wid∂Bϵ0=ϵ. Following [7], the coalescence flow is driven by the time evolution of ∂Bϵ0xs. One can therefore postulate that(13)Ordey·uxl,t~Ordey·uxr,tatxl,xr∈∂Bϵ0;t∈0,T.The placement of the velocity vector fields uxr,t, uxl,t is illustrated in Figure 3. The coalescence dynamics being dominated radially, consequently, the radial component of the equation of motion describing the evolution of the region Bϵ0xs—from the standpoint of Eulerian description—might at some time t>0 be written in the form(14)κ1∂t2b+κ2b2∂tb+κ3∂bγb=0fort∈0,T.The establishment of (14) is based on these chains of reasoning. In earlier sections, the standard basis ei:i=x,y,z was used. Here, use of the cylindrical basis eα:α=b,θ,z will be made to prove (14). Thus, every point x=x,y,zT∈Bϵ0xs on the basis ei:i=x,y,z can be represented by x=b,θ,yT∈Bϵ0xs on the basis eα:α=b,θ,z. For t∈0,T, it is postulated that the rate of change of linear momentum equates the differences of capillary and viscous forces across boundary ∂Bϵ0x0,i.e.,(15)RateofChangeofMomentum=∑α=cap,viscfα.With this view in mind, one then reasons as follows: The curve traced by a x,z-plane Πxzxs cutting the boundary ∂Bϵ0x0 is described by an equation of the form(16)zx=zs+b2-x-xs21/2on∂Bϵ0∩Πxz.Consequently, the curvature of that curve yields(17)curv∂Bϵ0x0=∂x2zx1+∂xzx23/2=1b.On the other hand, an ϵ-volume of Bϵ0x0 and an ϵ-surface of ∂Bϵ0x0 gives:(18)volϵBϵ0x0=ϵ·b·δb·δθandsurfϵ∂Bϵ0x0=ϵ·b·δθ,respectively. Note in passing that the plane Πxzxs and Πxzxs are clearly illustrated in Figure 4. Bearing in mind that the property b is a π;t-dependent scalar field, the rate of change of the linear momentum of volϵBϵ0 along eb, ϖBϵ0 say, then results in(19)ϖBϵ0=ρ∂t2b·volϵBϵ0eb.The capillary force difference sustained across the interface ∂Bϵ0x0 is given by(20)δfcap∂Bϵ0=-δp·surfϵ∂Bϵ0eb=-∂b2γ·curv∂Bϵ0·volϵBϵ0eb=-∂b2γb·volϵBϵ0eb.That due to viscous force is given by(21)δfvisc∂Bϵ0=∑i=b,θ,y-Sbiu·surfϵ∂Bϵ0b+Sbiu·surfϵ∂Bϵ0b+δbeb.=∑i=b,θ,y∂bbSbiu·volϵBϵ0ebIn (21), the tensor Su=(Siju)i,j=b,θ,y, with u-dependent entries Siju, where i is the row index and j the column index, stands for the viscous stress tensor. Substituting (19) through (21) into (15) yields(22)ρ∂t2b+∂b2γb-∑i=b,θ,y∂bbSbiu=0.For points x=b,θ,yT∈∂Bϵ0x0, it is asserted that u=ub,uθ,uyT~∂tb1,ϵ,ϵT on the basis of postulate (13), yielding(23)∑i=b,θ,y∂bbSbiu=μ1b∂b∂bub+1b2∂θ2ub+∂y2ub-ubb2-2b2∂θuθ~-μb2∂tb,sinceu=ur,uθ,uyT~∂tb1,ϵ,ϵT.Approximating the last term of (22) by (23) and then setting the triple κ1,κ2,κ3=ρ,μ,2, (14) follows. It is to be noted here that κii=13 is a π-dependent sequence because of the dependence of κ1 on ρ, and κ2 on μ, as is easily understood. One is now well prepared to discuss the so-called α-exponent laws, α reals.Figure 3
A freehanded sketch of the elementary domainBϵ0xs in R3.Figure 4
Illustration of the planesΠαzxs, the corresponding boundary curves Γαz (with α=x, y), and the location of the points xrear, xs, xfront∈⋃i=1,2‍Ωi∩S.Grantedbπ;t=ψπ·tα, after its substitution into (14), one arrives at the following algebraic expression:(24)αα-1κ1ψπtα-2+ακ2ψπt-α-1-κ3γψπ2t-2α=0.For that matter of analytical inquiries, set(25)Φ1π;t=αα-1κ1ψπtα-2,Φ2π;t=ακ2ψπt-α-1,Φ3π;t=-κ3γψπ2t-2αforgivenπ,andκii=13,t∈0,T,respectively. Mathematically, it is proposed to predict the power law behaviour of the property bπ;t based upon the magnitudes of Φii=13. Since, (24) is a three-term equation, the predictions are argued on three possible pairwise terms.


### 5.1. Case 1

Suppose thatOrdΦ1~OrdΦ2, and OrdΦ3≪1. The latter may then be neglected. Thus, for every fixed set of property constants, the following conditions must be fulfilled, namely, α-2=-α-1. The law describing this situation therefore writes(26)bπ;t=ψπ·t1/2fort∈0,T.


### 5.2. Case 2

Let it be assumed thatOrdΦ1~OrdΦ3, and OrdΦ2≪1, respectively. Then, analogous to the above statements, one obtains α-2=-2α. Thus, bπ;t is in this case of the form(27)bπ;t=ψπ·t2/3fort∈0,T.


### 5.3. Case 3

Conceive the orders of magnitudes ofΦii=13 with respect to themselves to be as follows: OrdΦ2~OrdΦ3, and OrdΦ1≪1. In this case, one deduces -α-1=-2α. Consequently, the scaling law describing the time evolution of the property bπ;t writes(28)bπ;t=ψπ·tfort∈0,T.It is to me pointed out that by OrdΦi is meant the order of magnitude of the property Φi (with i=1, 2, 3) with respect to time.For liquid droplets whose viscosities are not high enough, and having contact anglesθ<90∘, practical experiments have shown that bπ;t∝t2/3 [18, 19]. Contrarily, for liquids whose viscosities are high enough and having contact angles θ≳90∘, many authors have established the law bπ;t∝t1/2 [18, 21, 27]. On the other hand, in purely viscous coalescence regime, the law bπ;t∝t has also been deduced experimentally [3].In the next section are proved the power law regimes numerically, thereby giving agreement to experimental data.


## 5.1. Case 1

Suppose thatOrdΦ1~OrdΦ2, and OrdΦ3≪1. The latter may then be neglected. Thus, for every fixed set of property constants, the following conditions must be fulfilled, namely, α-2=-α-1. The law describing this situation therefore writes(26)bπ;t=ψπ·t1/2fort∈0,T.


## 5.2. Case 2

Let it be assumed thatOrdΦ1~OrdΦ3, and OrdΦ2≪1, respectively. Then, analogous to the above statements, one obtains α-2=-2α. Thus, bπ;t is in this case of the form(27)bπ;t=ψπ·t2/3fort∈0,T.


## 5.3. Case 3

Conceive the orders of magnitudes ofΦii=13 with respect to themselves to be as follows: OrdΦ2~OrdΦ3, and OrdΦ1≪1. In this case, one deduces -α-1=-2α. Consequently, the scaling law describing the time evolution of the property bπ;t writes(28)bπ;t=ψπ·tfort∈0,T.It is to me pointed out that by OrdΦi is meant the order of magnitude of the property Φi (with i=1, 2, 3) with respect to time.For liquid droplets whose viscosities are not high enough, and having contact anglesθ<90∘, practical experiments have shown that bπ;t∝t2/3 [18, 19]. Contrarily, for liquids whose viscosities are high enough and having contact angles θ≳90∘, many authors have established the law bπ;t∝t1/2 [18, 21, 27]. On the other hand, in purely viscous coalescence regime, the law bπ;t∝t has also been deduced experimentally [3].In the next section are proved the power law regimes numerically, thereby giving agreement to experimental data.


## 6. The Laws ofbπ;t Proved Numerically

Having established the power lawsbπ;t∝t1/2, t2/3, t theoretically, the chief aim of the present section is to validate them numerically. Before doing so, the groundwork necessary for such endeavours is first laid.Given the tripleπ-dependent proportionality constants ψαπα∈Iα, where Iα={1/2,2/3,1}. Write bπ;t=ψαπ·tα to stress that the prefactor referred to is the one associated with the α-exponent law and, let quantities btheo, bnum stand for its theoretical and numerical predictions, respectively. Thus, to express an ϵ-order discrepancy between the theoretical predictions and the simulation results, one can write btheo-bnum~Ordϵ and to negate the latter one can write btheo-bnum≠Ordϵ. The question arises: Which of those α-exponent laws, α∈Iα, should be computed first? The answer to this question lies in the following paragraph.Suppose two dropletsΩ1, Ω2 having the same size coalesce and time evolves, the bidroplet system ⋃i=1,2‍Ωi then will experience aviscous-inertial coalescence regime which is preceded by aviscous coalescence regime and followed by aninertial coalescence regime, as is easily understood. This shows the existence of two overlapping regimes. Moreover, if the set (0,Tα) designates the interval of time over which the law bπ;t=ψαπ·tα holds for some Tα>0, α∈Iα, then it is clear that the time intervals corresponding to those three coalescence regimes interrelate themselves as follows:(29)0,T1⊂0,T2/3⊂0,T1/2⊆0,T.These open sets are termed thecoalescence times of the corresponding α-regimes, α∈Iα. For some δTα>0, α∈Iα, it may be said that those two overlapping regimes occur over the time intervals(30)T1-δT1,T1+δT2/3andT2/3-δT1,T2/3+δT1/2.In view of the above discussion, it follows, in chronological order, that the establishment of 2/3-exponent law must precede that of the 1-exponent law and follow that of the 1/2-exponent law. Indeed, if one is to express for every α∈Iα the law bπ;t=ψαπ·tα by making use of the open set (0,Tα), the results would resolve into three algebraic expressions of the forms:(31)bπ;t≈b0+ψ1π·tfort∈0,T1,bπ;t≈b0+ψ1π·T1+ψ2/3π·t2/3fort∈T1,T2/3,bπ;t≈b0+ψ1π·T1+ψ2/3π·T2/32/3+ψ1/2π·t1/2fort∈T2/3,T1/2.Carrying out the discussion a step further, one can generalise the form of the scalar field bπ;t with respect to π;t. Supposed granted π=πii=1d, where d>1, and {ψαπ}α∈Iα a carefully chosen set of π-dependent proportionality constants, with the further supposition that Iα∈Q,i.e., Iα-{1/2,2/3,1}≠∅. Then, it is no error to state that the time evolution of bπ;t in the most general case is a series function of the form:(32)bπ;t≈b0+∑0<α≤1ψαπ·tαfort∈0,T.Leaving (32) out of consideration, thought it appears to demonstrate an important law, hence, might be a fairer, if not the fairest, approximation in practice, the following sections are devoted to the establishment of the laws bπ;t=ψαπ·tα, where α∈Iα.


### 6.1. The Viscous Regime Lawbπ;t=ψ1π·t

In this section, the numerical proof of the viscous coalescence power lawbπ;t=ψ1π·t is proved. How the numerical study has been undertaken may be understood through the following lines.Numerical experiments is first run. The geometrical properties used to construct the computational domain⋃i=1,2‍Ωi are those given in Table 1. The physical properties of the working fluid occupying the domain ⋃i=1,2‍Ωi as well as the initial and boundary conditions imposed on that domain are those tabulated in Table 1. The said power law is then established in the following manner:(33)bi=bπ;ti,ti=t0+i·δt,whereδt>0,i≥0,so that the pair of scalar fields ti,bii≥0 generates a sequence of points in the t,b-plane. Obviously, the dimensional units are such that [b]=[mm] and [t]=[s] (and [b]=[L·T-α], α∈Iα). For every i∈0,12, δt=5 [ms], the numerical value of ti,bi is extracted. This gives the necessary data required to describe the time evolution of the geometrical property bπ;t=ψ1π·t. By means of (33), the pairs of scalar fields ti,bii≥0 are collected in the viscous regime. When plotted in the t,b-plane, one obtains the graphical result exposed in Figure 5.Figure 5
Time evolution of the neck heightsbπ;t in the viscous regime. Thecircular markers∘ indicate the numerical experiment results and the solid line represents the fitted curve.In Figure5, circular markers ∘ are used to represent the numerical data, and the solid line the corresponding fitted power law. This power law function is shaped approximately as follows:(34)bπ;t=ψ1π·t≈0.0087t-4.9×10-6fort∈0,T.From Figure 5, one sees that btheo-bnum~Ordϵ only for t∈0,T1. The critical instant of time at which btheo-bnum≠Ordϵ is observed at t≈0.05 [s]; thus, approximately T1≈0.05 [s].


### 6.2. The Viscous-Inertial Regime Lawbπ;t=ψ2/3π·t2/3

Theviscous-inertial coalescence law bπ;t=ψ2/3π·t2/3 is demonstrated numerically in the present section. The procedure is similar to that discussed in the above section, differing only in the upper bound of the time T. That is, to demonstrate numerically the law bπ;t=ψ2/3π·t2/3, the property values exposed in Table 1 are used, and set T≈0.07 [s] (i.e., T>T1). Furthermore, since knowing the pair of scalar fields ti,bii≥0 yields the desired law in the t,b-plane numerically, to settle the 2/3-exponent law one then proceed in likely the same way as in the earlier section (see Section 6.1).In Figure6 are graphed the graphical results. Circular markers ∘ are used to represent the numerical data, and the solid line the equation corresponding fitted power law. That fitted power law curve is shaped as(35)bπ;t=ψ2/3π·t2/3≈0.0045t2/3fort∈0,T.It is understood that btheo-bnum~Ordϵ only for t∈0,T2/3. The critical instant of time at which btheo-bnum≠Ordϵ occurs at t≈0.06 [s]; thus, approximately T2/3≈0.06 [s].Figure 6
Time evolution of the neck heightsbπ;t in the inertial-viscous regime. Thecircular markers∘ indicate the numerical experiment results and the solid line represents the fitted curve.


### 6.3. The Inertial Regime Lawbπ;t=ψ1/2π·t1/2

The power lawbπ;t=ψ1/2π·t1/2 is here validated numerically. To do so numerical experiments were first run. The geometrical properties, initial and boundary conditions imposed, and the physical properties of the working fluid are those tabulated in Table 1. The said power law is then established in exactly the same way as in the foregoing two sections (see Sections 6.1 and 6.2).In Figure7 are exposed to view their graphical results. Circular markers ∘ are used to represent the numerical data, and the solid line the corresponding fitted power law. Through numerical experiments, it is believed that the following scalar function—representative of that solid line curve—gives good agreement:(36)bπ;t=ψ1/2π·t1/2≈0.0026t0.49fort∈0,T.In regard to (36), it follows that ψ1/2π≈0.0026. The power-the law bπ;t=ψ1/2π·t1/2 is thus proved numerically, demonstrated experimentally by [18, 21, 27]. It may be remarked that btheo-bnum~Ordϵ only for t∈0,T1/2. The critical instant of time at which btheo-bnum≠Ordϵ occurs at t≈0.17 [s]; thus, approximately T1/2≈0.17 [s]. The values of the scalar fields Ti>0 (i=1/2, 2/3, 1) are observed to be roughly equal to these values:(37)T1≈0.05s,T2/3≈0.06s,andT1/2≈0.17s.Moreover, from (36), (34), (35), one also observes that(38)
O
r
d
ψ
α
π
~
1
0
-
3
m
m
·
s
-
α
f
o
r
a
l
l
α
∈
I
α
.In what follows the π-dependent sequence ψαπα∈Iα are analysed.Figure 7
Time evolution of the neck heightsbπ;t in the inertial regime. Thecircular markers∘ indicate the numerical experiment results and the solid line represents the fitted curves.


### 6.4. Of the Characters ofψαπα∈Iα

When the literature was surveyed, one found no more than two elements of the sequenceψαπα∈Iα have gained attention, yet their establishment in terms of the property constants involved in the problems they ultimately lie have been correctly settled. With this view in mind, it is proposed to discuss in the present section the characters of ψαπα∈Iα through analysis.The analysis commences by relating the leading property associated with the interfacial coalescence hydrodynamics to that associated with the bulk coalescence hydrodynamics, that is, it is proposed to relatebπ;t to ux,t on Γxz. It is clear that bπ;t=x-xs, where x, xs∈Πxzxs. Therefore, taking its derivative, one obtains(39)u·x-xs=αψαπ2t2α-1onΓxz×0,Tα,whereα∈Iα.On the other hand, the leading property constants being ρ, μ, γ0, b0, R1, consequently, π=ρ,μ,γ0,b0,R1T. Given the characteristic scales Tα,Xα,Uα,Ψαα∈Iα, one rescales the properties appearing in (39) as follows: t↦Tαt, x↦xs+Xαx-xs, u↦Uαu, and ψα↦Ψαψα, yielding(40)OrdΨαπ~UαXααTα2α-1,whereα∈Iα.It now remains to express the right-hand side of (40) in terms of the property vector π=ρ,μ,γ0,b0,R1T. To do this, one reasons according to the α-regimes, where α∈Iα. Analysing the set of properties Tα,Xα,Uαα∈Iα in terms of π, the following orders of magnitude follow:(41)OrdT1,T2/3,T1/2~μb0γ0,μb0γ0,ρR13γ0,OrdX1,X2/3,X1/2~b0,b0,R1,OrdU1,U2/3,U1/2~γ0μ,γ0μ,γ0ρR1.These being so, substituting (41) into (40), one finds(42)OrdΨ1π,Ψ2/3π,Ψ1/2π~γ0μ,32γ02b0μ21/3,2γ0R1ρ1/4.Moreover, inserting the property constants of Table 1 into (42) and, thereafter, dividing each element of (ψαπtheo)α∈Iα by the corresponding element of (ψαπnum)α∈Iα, the orders of magnitude of the ratio of Ψαπtheo to Ψαπnum yield(43)OrdΨαπtheoΨαπnum≈4.5~Ord1forallα∈Iα.In (43), it is understood that quantities Ψαπtheo, Ψαπnum stand for the theoretical and numerical estimation, respectively, of the scalar field Ψαπ. Expressing (34), (35) and (36) in terms of (42) gives(44)bπ;t~γ0μtfort∈0,T1,bπ;t~32γ02b0μ21/3t2/3fort∈0,T2/3,bπ;t~2γ0R1ρ1/4t1/2fort∈0,T1/2.Equation (44) at once justifies the reasonableness of the foregoing chains of reasoning and, this completes the discussion of the present section concerning the π;t-dependent evolution of the field bπ;t. In the next section, it is proposed to enter upon the discussion of the characters of the curve ux,t along a carefully chosen one-dimensional boundary.


### 6.5. Of the Characters ofux,t

If the velocity normu of the interacting droplets systems is plotted on the curved boundary Γyz, the curve traced by the intersection of the interface ∂Ω1∪Ω2 with the plane Πyzxs, one obtains Figure 8. The placement of arrows on the u curves of Figure 8 serve merely to indicate the evolution of time evolution of these curves. In Figure 8, the graphs of the velocity norm u correspond to the sequence of times tkk=05=(10k) [ms]. The nature of these curves are now explained.Figure 8
Time evolution of the velocity normu of two interacting droplets with radii R1, R2 (with R2=R1) along the boundary curve Γyz for the sequence of times tkk=05=(10k)k=05 [ms].Inspecting Figure8, one sees that the highly curved meniscus induces an ϵ-order singularity into the coalescence flow in the early-time regime, peaking to magnitudes of the order of u≈20 mm·s-1 inside that highly curved meniscus zone (indicated by dashed, vertical lines). It is perhaps because of this striking effect that some authors have characterised such motion assingular motion [7, 22].It is interesting to see the spatiotemporal distributions of the field variableux,t. This is the purpose of the next section.


## 6.1. The Viscous Regime Lawbπ;t=ψ1π·t

In this section, the numerical proof of the viscous coalescence power lawbπ;t=ψ1π·t is proved. How the numerical study has been undertaken may be understood through the following lines.Numerical experiments is first run. The geometrical properties used to construct the computational domain⋃i=1,2‍Ωi are those given in Table 1. The physical properties of the working fluid occupying the domain ⋃i=1,2‍Ωi as well as the initial and boundary conditions imposed on that domain are those tabulated in Table 1. The said power law is then established in the following manner:(33)bi=bπ;ti,ti=t0+i·δt,whereδt>0,i≥0,so that the pair of scalar fields ti,bii≥0 generates a sequence of points in the t,b-plane. Obviously, the dimensional units are such that [b]=[mm] and [t]=[s] (and [b]=[L·T-α], α∈Iα). For every i∈0,12, δt=5 [ms], the numerical value of ti,bi is extracted. This gives the necessary data required to describe the time evolution of the geometrical property bπ;t=ψ1π·t. By means of (33), the pairs of scalar fields ti,bii≥0 are collected in the viscous regime. When plotted in the t,b-plane, one obtains the graphical result exposed in Figure 5.Figure 5
Time evolution of the neck heightsbπ;t in the viscous regime. Thecircular markers∘ indicate the numerical experiment results and the solid line represents the fitted curve.In Figure5, circular markers ∘ are used to represent the numerical data, and the solid line the corresponding fitted power law. This power law function is shaped approximately as follows:(34)bπ;t=ψ1π·t≈0.0087t-4.9×10-6fort∈0,T.From Figure 5, one sees that btheo-bnum~Ordϵ only for t∈0,T1. The critical instant of time at which btheo-bnum≠Ordϵ is observed at t≈0.05 [s]; thus, approximately T1≈0.05 [s].


## 6.2. The Viscous-Inertial Regime Lawbπ;t=ψ2/3π·t2/3

Theviscous-inertial coalescence law bπ;t=ψ2/3π·t2/3 is demonstrated numerically in the present section. The procedure is similar to that discussed in the above section, differing only in the upper bound of the time T. That is, to demonstrate numerically the law bπ;t=ψ2/3π·t2/3, the property values exposed in Table 1 are used, and set T≈0.07 [s] (i.e., T>T1). Furthermore, since knowing the pair of scalar fields ti,bii≥0 yields the desired law in the t,b-plane numerically, to settle the 2/3-exponent law one then proceed in likely the same way as in the earlier section (see Section 6.1).In Figure6 are graphed the graphical results. Circular markers ∘ are used to represent the numerical data, and the solid line the equation corresponding fitted power law. That fitted power law curve is shaped as(35)bπ;t=ψ2/3π·t2/3≈0.0045t2/3fort∈0,T.It is understood that btheo-bnum~Ordϵ only for t∈0,T2/3. The critical instant of time at which btheo-bnum≠Ordϵ occurs at t≈0.06 [s]; thus, approximately T2/3≈0.06 [s].Figure 6
Time evolution of the neck heightsbπ;t in the inertial-viscous regime. Thecircular markers∘ indicate the numerical experiment results and the solid line represents the fitted curve.


## 6.3. The Inertial Regime Lawbπ;t=ψ1/2π·t1/2

The power lawbπ;t=ψ1/2π·t1/2 is here validated numerically. To do so numerical experiments were first run. The geometrical properties, initial and boundary conditions imposed, and the physical properties of the working fluid are those tabulated in Table 1. The said power law is then established in exactly the same way as in the foregoing two sections (see Sections 6.1 and 6.2).In Figure7 are exposed to view their graphical results. Circular markers ∘ are used to represent the numerical data, and the solid line the corresponding fitted power law. Through numerical experiments, it is believed that the following scalar function—representative of that solid line curve—gives good agreement:(36)bπ;t=ψ1/2π·t1/2≈0.0026t0.49fort∈0,T.In regard to (36), it follows that ψ1/2π≈0.0026. The power-the law bπ;t=ψ1/2π·t1/2 is thus proved numerically, demonstrated experimentally by [18, 21, 27]. It may be remarked that btheo-bnum~Ordϵ only for t∈0,T1/2. The critical instant of time at which btheo-bnum≠Ordϵ occurs at t≈0.17 [s]; thus, approximately T1/2≈0.17 [s]. The values of the scalar fields Ti>0 (i=1/2, 2/3, 1) are observed to be roughly equal to these values:(37)T1≈0.05s,T2/3≈0.06s,andT1/2≈0.17s.Moreover, from (36), (34), (35), one also observes that(38)
O
r
d
ψ
α
π
~
1
0
-
3
m
m
·
s
-
α
f
o
r
a
l
l
α
∈
I
α
.In what follows the π-dependent sequence ψαπα∈Iα are analysed.Figure 7
Time evolution of the neck heightsbπ;t in the inertial regime. Thecircular markers∘ indicate the numerical experiment results and the solid line represents the fitted curves.


## 6.4. Of the Characters ofψαπα∈Iα

When the literature was surveyed, one found no more than two elements of the sequenceψαπα∈Iα have gained attention, yet their establishment in terms of the property constants involved in the problems they ultimately lie have been correctly settled. With this view in mind, it is proposed to discuss in the present section the characters of ψαπα∈Iα through analysis.The analysis commences by relating the leading property associated with the interfacial coalescence hydrodynamics to that associated with the bulk coalescence hydrodynamics, that is, it is proposed to relatebπ;t to ux,t on Γxz. It is clear that bπ;t=x-xs, where x, xs∈Πxzxs. Therefore, taking its derivative, one obtains(39)u·x-xs=αψαπ2t2α-1onΓxz×0,Tα,whereα∈Iα.On the other hand, the leading property constants being ρ, μ, γ0, b0, R1, consequently, π=ρ,μ,γ0,b0,R1T. Given the characteristic scales Tα,Xα,Uα,Ψαα∈Iα, one rescales the properties appearing in (39) as follows: t↦Tαt, x↦xs+Xαx-xs, u↦Uαu, and ψα↦Ψαψα, yielding(40)OrdΨαπ~UαXααTα2α-1,whereα∈Iα.It now remains to express the right-hand side of (40) in terms of the property vector π=ρ,μ,γ0,b0,R1T. To do this, one reasons according to the α-regimes, where α∈Iα. Analysing the set of properties Tα,Xα,Uαα∈Iα in terms of π, the following orders of magnitude follow:(41)OrdT1,T2/3,T1/2~μb0γ0,μb0γ0,ρR13γ0,OrdX1,X2/3,X1/2~b0,b0,R1,OrdU1,U2/3,U1/2~γ0μ,γ0μ,γ0ρR1.These being so, substituting (41) into (40), one finds(42)OrdΨ1π,Ψ2/3π,Ψ1/2π~γ0μ,32γ02b0μ21/3,2γ0R1ρ1/4.Moreover, inserting the property constants of Table 1 into (42) and, thereafter, dividing each element of (ψαπtheo)α∈Iα by the corresponding element of (ψαπnum)α∈Iα, the orders of magnitude of the ratio of Ψαπtheo to Ψαπnum yield(43)OrdΨαπtheoΨαπnum≈4.5~Ord1forallα∈Iα.In (43), it is understood that quantities Ψαπtheo, Ψαπnum stand for the theoretical and numerical estimation, respectively, of the scalar field Ψαπ. Expressing (34), (35) and (36) in terms of (42) gives(44)bπ;t~γ0μtfort∈0,T1,bπ;t~32γ02b0μ21/3t2/3fort∈0,T2/3,bπ;t~2γ0R1ρ1/4t1/2fort∈0,T1/2.Equation (44) at once justifies the reasonableness of the foregoing chains of reasoning and, this completes the discussion of the present section concerning the π;t-dependent evolution of the field bπ;t. In the next section, it is proposed to enter upon the discussion of the characters of the curve ux,t along a carefully chosen one-dimensional boundary.


## 6.5. Of the Characters ofux,t

If the velocity normu of the interacting droplets systems is plotted on the curved boundary Γyz, the curve traced by the intersection of the interface ∂Ω1∪Ω2 with the plane Πyzxs, one obtains Figure 8. The placement of arrows on the u curves of Figure 8 serve merely to indicate the evolution of time evolution of these curves. In Figure 8, the graphs of the velocity norm u correspond to the sequence of times tkk=05=(10k) [ms]. The nature of these curves are now explained.Figure 8
Time evolution of the velocity normu of two interacting droplets with radii R1, R2 (with R2=R1) along the boundary curve Γyz for the sequence of times tkk=05=(10k)k=05 [ms].Inspecting Figure8, one sees that the highly curved meniscus induces an ϵ-order singularity into the coalescence flow in the early-time regime, peaking to magnitudes of the order of u≈20 mm·s-1 inside that highly curved meniscus zone (indicated by dashed, vertical lines). It is perhaps because of this striking effect that some authors have characterised such motion assingular motion [7, 22].It is interesting to see the spatiotemporal distributions of the field variableux,t. This is the purpose of the next section.


## 7. The Coalescence Flow Fields Illustrated

Having discussed the growth of the time-dependent functionbπ;t in the inertial, inertial-viscous and viscous regimes, the evolutions of the vector field function ux,t with respect to x,t∈⋃i=1,2‍Ωi×0,T are now reported. In few words, the present endeavours is an attempt to lay additional stress on the qualitative aspects of the proposed model.When the coalescence domain⋃i=1,2‍Ωi undergoes coalescence phenomenon, as soon as its pair of droplets are brought into contact, a mass movement is started, thereby causing the corresponding droplets to fuse. It is not necessary to go into details; sufficient is to demonstrate that the proposed model is capable of describing the underlying motion phenomena and the latter are exposed in the 3D space and on the 2D planes Παzxs (with α=x, y).In Figures9 and 10, the time evolutions of the velocity norm u of the coalescence domain ⋃i=1,2‍Ωi at times t=5, 50, 100, 190 [ms] are exposed in the 3D space. The corresponding two-dimensional plots of the scalar field u are reported, respectively, on the 2D planes Παzxs (with α=x, y) given in Figures 11–14, where solid lines stand for the iso-u contours.Figure 9
Three-dimensional plots of the velocity normu of two interacting droplets of the same radius R1=R2 at times t=5, 50 [ms].Figure 10
Three-dimensional plots of the velocity normu of two interacting droplets of the same radius R1=R2 at times t=100, 190 [ms].Figure 11
Two-dimensional plots of the velocity normu of two interacting droplets of the same radius R1=R2 on the plane Πxzxs at times t=5, 50 [ms]. The solid lines stand for the iso-u contours.Figure 12
Two-dimensional plots of the velocity normu of two interacting droplets of the same radius R1=R2 on the plane Πxzxs at times t=100, 190 [ms]. The solid lines stand for the iso-u contours.Figure 13
Two-dimensional plots of the velocity normu of two interacting droplets of the same radius R1=R2 on the plane Πxzxs at times t=5, 50 [ms]. The solid lines stand for the iso-u contours.Figure 14
Two-dimensional plots of the velocity normu of two interacting droplets of the same radius R1=R2 on the plane Πxzxs at times t=100, 190 [ms]. The solid lines stand for the iso-u contours.


## 8. Closing Remarks

In this part, a three-dimensional modelling methodology has been proposed to study droplets coalescence phenomenon. On the whole, the statements arrived at conclusively are:When, upon variations of geometrical parameters, the footprint⋃i=1,2‍Ωi∩S and the free-surface ∂(Ω1∪Ω2) of the bidroplet system ⋃i=1,2‍Ωi reduce to the footprint Ω2∩S and the free-surface ∂Ω2, respectively, of a monodroplet system Ω2, the corresponding boundary conditions are by default adapted into the COMSOL interfaces. As a result, with one numerical model, its property constants can be so arranged so as to compute many physical situations, a task which one would have undertaken had the proposed model been unversatile.Taking advantage of the geometrical features of the proposed model, it is possible to deduce through theoretical analysis the scaling lawsbπ;t=ψαπ·tα (with α=1/2, 2/3, 1), laws already reported in the (experimental) literature.Carried out numerical proofs of these laws,α≈1 is observed in the visco-dominated coalescence flow, α≈2/3 in the inertia-visco-dominated coalescence flow, and α≈1/2 in the inertia-dominated coalescence flow.Through scaling analysis, it is found that the dependence of the sequence of prefactors(ψαπ)α∈Iα on π=ρ,μ,γ0,b0,R1T can be approximated as follows: ψ1π≈γ0/μ, ψ2/3π≈3/2(γ02b0/μ2)1/3 and ψ1/2π≈2(γ0R1/ρ)1/4.The flow field variables illustrated by2D, 3D plots clearly demonstrate that the expected features associated with droplets coalescence phenomenon are well captured by the model.On the whole, yet the approach in the present analysis has been to present the leading features of the proposed model by considering the coalescence that results under an idealized situation, it is no way contradictory to declare that it can be used to model more complex physical situations. For example, numerical experiments on droplets coalescence phenomenon can be carried out under appropriate CFD conditions to display useful information concerning mixing process driven by Marangoni flow.



---

*Source: 4906016-2018-10-18.xml*