doi:10.1017/jfm.2016.381
Controlling droplet bouncing and coalescence with surfactant
K.-L. Pan1,†, Y.-H. Tseng2, J.-C. Chen1, K.-L. Huang1, C.-H. Wang1 and M.-C. Lai3
1Department of Mechanical Engineering, National Taiwan University, Taipei 10617, Taiwan, ROC 2Department of Applied Mathematics, National University of Kaohsiung, Kaohsiung 81148,
Taiwan, ROC
3Department of Applied Mathematics, National Chiao Tung University, Hsinchu 30050, Taiwan, ROC (Received 6 April 2015; revised 28 March 2016; accepted 31 May 2016)
The collision between aqueous drops in air typically leads to coalescence after impact. Rebounding of the droplets with similar sizes at atmospheric conditions is not generated, unless with significantly large pressure or high impact parameters exhibiting near-grazing collision. Here we demonstrate experimentally the creation of a non-coalescent regime through addition of a small amount of water-soluble surfactant.
We perform a direct simulation to account for the continuum and short-range flow dynamics of the approaching interfaces, as affected by the soluble surfactant. Based on the immersed-boundary formulation, a conservative scheme is developed for solving the coupled surface-bulk convection–diffusion concentration equations, which presents excellent mass preservation in the solvent as well as conservation of total surfactant mass. We show that the Marangoni effect, caused by non-uniform distributions of surfactant on the droplet surface and surface tension, induces stresses that oppose the draining of gas in the interstitial gap, and hence prohibits merging of the interfaces.
In such gas–liquid systems, the repulsion caused by the addition of surfactant, as frequently observed in liquid–liquid systems such as emulsions in the form of an electric double-layer force, was found to be too weak to dominate in the attainable range of interfacial separation distances. These results thus identify the key mechanisms governing the impact dynamics of surfactant-coated droplets in air and imply the potential of using a small amount of surfactant to manipulate impact outcomes, for example, to prevent coalescence between droplets or interfaces in gases.
Key words: breakup/coalescence, drops, gas/liquid flows
1. Introduction
The collision dynamics between two droplets in the gas phase plays a crucial role in various disciplines of nature and practical interest, such as the formation of raindrops (Gunn 1965), the operation of nuclear reactors (Bauer, Bertsch & Schulz 1992; Moretto et al. 1992), spray combustion (Chiu 2000) and fire-fighting with
† Email address for correspondence: [email protected]
(I) Coalescence B
(III) Coalescence
(V) Stretching separation
(IV) Reflexive separation (II) Bouncing
FIGURE 1. Typical diagram showing the regimes of collision outcomes (Qian &
Law 1997). I: coalescence with minor deformation; II: bouncing; III: coalescence with substantial deformation; IV: near head-on, also known as reflexive separation (Ashgriz &
Poo1990); V: off-centred, also known as stretching separation. Here WeS and WeH indicate the soft and hard transition boundary, respectively, and Wec designates the critical value beyond which the merged droplets separate after coalescence of two droplets with head-on impact.
liquid injection (Grant, Brenton & Drysdale 2000). Under the complication due to a large number of key factors leading to redistribution of liquid mass, momentum and energy, an essential issue is how to manipulate the impact of droplets so as to achieve the desired performance. This mostly concerns whether coalescence occurs. As a rudimentary mechanism for droplet collision, the generic scenarios after the impact of two droplets made of identical liquid and size in gases have been widely investigated (Jiang, Umemura & Law 1992; Qian & Law 1997; Pan, Law & Zhou 2008). They show characteristic transitions from (I) coalescence after minor deformation to (II) bouncing to (III) coalescence after substantial deformation and to (IV) temporary coalescence, followed by the separation of primary and satellite droplets. Further breakup and splattering of the merged droplet into multiple secondary droplets can be created when the impact energy is relatively large (Pan, Chou & Tseng 2009), similar to the splashing phenomena observed in droplet impact upon a surface (Yarin 2006). As shown in figure 1 (Qian & Law 1997), such a scenario is generated when the key parameter, the Weber number, is increased, as defined by We = ρiU∞2D0/σ0, which indicates the ratio between the kinetic energy and the surface energy, where U∞ is the relative impact speed, D0 is the diameter of the droplet and ρi and σ0, respectively, are the density and surface tension of the liquid. Here WeS indicates the soft transition boundary when the impact is characterized by minor deformation, and WeH is the hard transition boundary where substantial deformation is created, as analysed in Pan et al. (2008). Wec designates the critical value (Ashgriz & Poo 1990; Jiang et al. 1992) beyond which the merged droplets separate after temporary coalescence of two droplets in head-on impact. The typical sequences of the processes are presented in figure 2, which are to be discussed further in §4 along with our experimental findings on the effects of adding surfactant.
It is generally known (Adam, Lindblad & Hendricks 1968; Brazier-Smith, Jennings
& Latham 1972; Ashgriz & Poo 1990) that bouncing is not created during collision under atmospheric conditions between two aqueous droplets of identical size, although it may happen if the environmental pressure is significantly larger (Qian & Law
0.0 ms 0.2 ms (a)
(c) (d)
(e)
(b)
0.4 ms
0.0 ms 0.2 ms 0.4 ms 0.6 ms
0.8 ms 1.4 ms 1.6 ms 2.2 ms
2.6 ms 3.0 ms 3.2 ms
0.0 ms 0.2 ms 0.4 ms
0.6 ms 1.2 ms 2.2 ms
2.6 ms 2.8 ms 3.0 ms
3.4 ms
0.0 ms 0.2 ms 0.4 ms 0.6 ms
0.8 ms 1.2 ms 1.4 ms 1.6 ms
2.0 ms 2.6 ms 3.0 ms 3.2 ms
0.6 ms 0.8 ms 1.0 ms
0.0 ms 0.2 ms 0.4 ms
0.6 ms 0.8 ms 1.0 ms
FIGURE 2. Characteristic sequences of droplet impact experiments: (a) coalescence with minor deformation, regime I (We = 0.72, B = 0.026); (b) coalescence with substantial deformation, regime III (We = 11.77, B = 0.050); (c) bouncing at large impact parameter, regime II (We = 13.10, B = 0.585); (d) stretching separation at large impact parameter, regime V (We = 25.32, B = 0.586); (e) reflexive separation at small impact parameter, regime IV (We = 26.20, B = 0.003).
1997). It can also occur when the impact parameter, B = χ/D0, is so large that the collision path deviates significantly from head-on (B = 0) and approaches the grazing condition (B = 1), where χ is the projection of the distance between the droplet centres in the direction normal to the velocity vector, U∞. A similar phenomenon of non-coalescence is also observed in the oblique collision of fluid jets (Wadhwa, Vlachos & Jung 2013). Phenomenologically, the propensity for bouncing or merging is a consequence of the ease with which the gas in the interdroplet gap can be squeezed out by the colliding interfaces such that they can make contact at the molecular level, leading to their destruction and consequent merging (Pan et al. 2008). The essential role of the interdroplet gaseous film has been experimentally demonstrated through the dependence of the collision outcome on the density of the gas medium (Qian
& Law 1997). This is analogous to the air cushioning effect identified recently for
single droplet impact upon a solid surface, as discussed in Hicks & Purvis (2010) and Kolinski et al. (2012).
To manipulate the collision consequence without substantial modification of the system, one could vary the fluid properties by additives. Without changing the bulk fluid properties much, controlling the transition between the rebounding and the deposition regimes of an aqueous droplet impacting a solid surface has been studied by adding a small amount of polymer, as reported in Bergeron et al. (2000). In our recent work (Pan & Hung 2010) regarding collision between a droplet and a wet surface, addition of a surfactant was found to modify ostensibly the transformation boundaries on the regime diagram. Surfactant means surface active agent and tends to attach to the surface rather than dissolve in the bulk fluid. Its effect on the fluid can be distinct from other additives that merely change the fluid properties, such as glycerol which varies viscosity, since the characteristics around the interfaces can be changed remarkably. Although such studies have not been performed for binary droplet impact in gases, a surfactant is frequently employed in liquid–liquid systems to prevent coalescence of droplets in liquids. The physical mechanisms are attributed to various factors, such as non-uniform distribution of the surfactant concentration that yields Marangoni stresses (Yeo et al. 2003; Dai & Leal 2008) and repulsive intermolecular forces (Petsev 2000), specifically the electric double-layer (EDL) force (Zhang et al. 2010), but their exact roles in the processes are not fully understood.
In the present liquid–gas configuration without electric charges in the interstitial film between the impact droplets, influence of such short-range forces can be essentially neglected, so the Marangoni effect can be tested for its role exclusively in the collision process. Through a systematic study using experiments on the impact between two water droplets of identical size in air, accounting for variations of We and B, and the full Navier–Stokes computations including the effects of soluble surfactant and intermolecular forces, we present for the first time the detailed flow dynamics and correlations of these factors. Our results indicate substantial enhancement of droplet rebounding with addition of surfactant and show unambiguously the dominance of the Marangoni effect in the process. By adding a small amount of surfactant, the collision outcomes could be manipulated in the desired manner, e.g. to enhance the stability of liquid–gas systems by preventing coalescence of droplets.
Solving numerically the present problems of droplet impact is challenging due to the need for high resolution at the impact region as well as the coupling of surfactant equations (convection–diffusion equations) with moving interfaces (whereby surfactant is contained only inside the droplets). The resolution issue can be handled by either the strategy of a fixed structure grid (Vinokur 1983) or adaptive mesh refinement (AMR) techniques (Berger & Oliger 1984; Berger & Colella 1989). While how to solve the surfactant equations efficiently and how to preserve the surfactant exactly in one of the phases i.e. conservation of surfactant mass, are the main difficulties for research interest, there have been achievements in such simulations, although mostly regarding insoluble surfactants. Some numerical methods such as the surface element method (Burger 2005; Dziuk & Elliott 2007), level set method (Bertalmio et al. 2001; Adalsteinsson & Sethian 2003; Leung, Lowengrub & Zhao 2011) and phase field method (Rätz & Voigt 2006; Teigen et al. 2009; Elliott et al. 2011), have been developed to solve the convection–diffusion equations with evolving interfaces.
The front tracking method (Peskin 1972, 2002; Unverdi & Tryggvason 1992) is generally more accurate for tracking interfaces, but more complicated implementation is needed for restructuring of surface mesh, specifically when it suffers significant deformation or even topological changes. As reported in Lai, Tseng & Huang (2008),
we have successfully developed a mass-conservative scheme for problems of moving interfaces with convection–diffusion equations and applied it to simulate interfacial flows with insoluble surfactant (Lai et al. 2008; Lai, Tseng & Huang 2010). Recently, a conservative numerical method for soluble surfactant cases (Chen & Lai 2014) was developed to study droplet deformation in a contaminated two-dimensional fluid system. Surfactant concentrations both in the bulk fluid and on the interface were considered, and a simple adsorption–desorption model was adopted to describe the exchange of surfactant at an interface. In the present work, we extend the methodology described in Chen & Lai (2014) for an axisymmetric system and investigate how the surfactant affects the head-on collision of two identical droplets in a gaseous environment. We formulate the coupled surface-bulk convection–diffusion equations in the framework of immersed boundaries (Peskin 1977) so that the adsorption and desorption processes can be termed as a singular source in the bulk equation.
Moreover, by using an indicator function, we can embed the bulk equation in the whole computational domain; thereby a regular Eulerian finite-difference scheme can be implemented straightforwardly without further concern for the complexity when treating the moving surfaces as irregular boundaries in the domain of computation.
The numerical scheme can preserve the total mass of surfactant exactly in a discrete sense. By introducing the indicator function and solving the bulk equation in the whole computational domain, one can avoid evaluating the surfactant flux across the interface due to adsorption and desorption processes.
To describe the work according to the findings in the experiments and the computational analyses used to interpret the mechanisms, the report is organized as follows. In the next section, the experimental approach is illustrated. Since similar methodologies have been described in our previous studies (Pan et al. 2008, 2009), not much exposition is intended. Due to the specific numerical simulation approach on the other hand, more details are reported thereafter for the computational analyses.
In §3, the mathematical model based on the Navier–Stokes equations and surfactant concentration equations is presented, including the corresponding interfacial forces and intermolecular forces. A strictly mass-conservative immersed-boundary method is described in the section. The general findings in experiments with addition of surfactant are then discussed in §4. This is followed by §5, in which a series of parametric studies via numerical simulation are performed to identify the significance of various competing physical mechanisms in the process of droplet rebounding.
Specifically, the dominance of tangential stresses is verified, and the features of surfactant solubility are discussed as well. Moreover, the influence of viscosity is investigated in the section, both experimentally and computationally. The conclusions are summarized in §6.
2. Experimental set-up
The generation of the droplets for the desired collision condition is based on the drop-on-demand method that is similar to ink printing technologies (Ashgriz
& Poo 1990; Jiang et al. 1992; Qian & Law 1997; Pan et al. 2008). As shown in figure 3, two droplets of identical size and material are generated by nozzles triggered by the vibration of piezoelectric plates. They are made to impinge onto each other in a controlled path with adjustable angles. Time-resolved images are either taken by a standard CCD through stroboscopy synchronized with the droplet generation circuit or recorded by a high-speed CMOS digital camera (X-StreamTM Vision, XS-4), which supports a maximum resolution of 512 × 512 pixels with
1 4
5
3
8 9
10 2
6
7
FIGURE 3. (Colour online) Schematic diagram of the experimental set-up with components 1: desktop, 2: LED-based illuminator, 3: high-speed digital camera, 4:
oscilloscope, 5: electronic control box, 6: enclosure of test section, 7: liquid supply reservoir, 8: motion controller for droplet generator, 9: droplet generator, 10: liquid supply tube.
5100 frames per second (fps). Depending on the test conditions, the frame rate can be raised to 20 000 f.p.s while the resolution is reduced to 128 × 512. The shutter of the high-speed camera is synchronized with an LED lamp that can support the shortest duration of 1 µs to capture images with a sufficiently small exposure and adequate light intensity while avoiding blurring due to background scattering. A digital imaging system then accurately time resolves the collision event, records the droplet image and processes the data. These approaches provide image capturing suitable for various conditions and can yield fine temporal and spatial resolutions for the transitional behaviours between different regimes. In addition, the boundaries of the droplets in the images are detected based on the Hough transform (Illingworth & Kittler 1988) coded by a Matlab program (as downloaded from http://www.mathworks.com/matlabcentral/fileexchange/9168, which was developed by Tao Peng in 2005 and updated in 2010). The Hough transform is a technique useful in image processing and computer vision, which can be used to isolate features of specific shape within an image. This method is advantageous for the identification of complex structure with relatively low sensitivity to image noise or missing data. Via the detection for the sharpest gradient of greyscale, the droplet shape and hence the diameter, moving distance and velocity can be measured with high accuracy. More details of the experimental methodology have been given in Chou (2008) and Chen (2010).
Three primary categories of surfactants have been tested in the experiments (Pan &
Chen 2012), i.e. anionic (S111n), amphoteric (S131) and non-ionic (S386), which are
0 10 20 30 40 50 60 70 80
S111n S131 S386
0.5 1.0 1.5
FIGURE 4. (Colour online) Variation of measured surface tension with surfactant concentration. The CMC limit of S111n is defined as ψ0∼ 0.35 %.
all soluble in water. The density and viscosity of the aqueous solutions are essentially identical, but the surface tension (σ0) changes significantly, as shown in figure 4, with the variation of initial surfactant concentration, ψ0. It is seen that σ0 attains a minimum at the critical micelle concentration (CMC) limit, which designates the saturation level of surfactant on the liquid interface. These surfactants belong to a series of Surflon surfactants, produced by AGC Seimi Chemical Co., LTD. The surfactant, with a purity of 30 %, is a fluorochemical surface acting agent composed of perfluoroalkyl betaine. This product possesses excellent attributes that are not possible for hydrocarbon-based surfactants, such as a significant reduction in the surface tension with a small quantity added. While its viscosity is 8.8 cP (10−3 Pa s), the addition of a small amount does not substantially change the viscosity of its water solution. The properties of fluids tested in the study are listed in table 1. It is seen that, while surface tension is reduced ostensibly, viscosity is almost invariant after the surfactant is added to water. In contrast, when glycerol is added, the viscosity increases prominently while the surface tension stays essentially the same and density increases slightly.
In passing we note that, due to more data available for S111n in the broader range of transition before reaching the CMC limit, and hence more precise control of surface tension with addition of surfactant, as seen in the measurement in figure 4, most of the results presented are based on the surfactant S111n. While the ions on the droplet interfaces might induce certain effects, such a factor is shown to be insignificant for the key phenomena to be discussed in the following.
3. Mathematical model and numerical method
Numerical simulations are performed for multiphase fluid dynamics based on the immersed-boundary (IB) method (Peskin 1977). With conservation of specific mathematical and physical quantities, the surface forces are smeared to surrounding fluid grid points, so the computation can be carried out in a unified domain. This is generally known as an Eulerian–Lagrangian approach (Kuan, Pan & Shyy 2014)
Surfactant weight (wt) Density Viscosity Surface tension Uncertainty concentration ψ0 ( %) ρi (kg m−3) µi (cP) σ0 (N m−1) 1σ0/σ0 ( %)
Water 0 998 1.02 0.072 1.25
S111n 0.1 998 1.02 0.032 3.12
S131 0.005 998 1.02 0.032 3.44
S386 0.005 998 1.02 0.032 3.53
S111n 0.3 998 1.02 0.021 5.98
S111n 1.0 998 1.02 0.017 5.59
Glycerol 0 1069 2.83 0.071 1.33
(wt 30 %)
S111n+ 1.0 1069 2.83 0.017 5.88
Glycerol (wt 30 %)
Ethanol 0 789 1.20 0.022 1.85
TABLE 1. Properties of the tested fluids.
for moving-boundary problems. In addition, the dynamics of soluble surfactant is calculated along with the conservation equations of fluid momentum and mass.
3.1. Governing equations
In this subsection, the mathematical model for head-on collision of two equal sized droplets with soluble surfactant in a three-dimensional (3-D) axisymmetric immersed-boundary formulation is described. We assume that, in the computational domain Ω = {(r, z)|r ∈ (0, Lr),z ∈ (0, Lz)}, the two droplets are placed symmetrically with respect to z = 0 (r axis) so that the collision direction is along the z axis. The fluid domain Ω consists of the droplet phase Ωi and the surrounding fluid Ωo (gas, throughout the report), and so in between there is an interface Σ separating them.
Mathematically, this interface is the boundary of Ωi, i.e. Σ = ∂Ωi. In this work, for insoluble surfactant, we mean that the surfactant exists only on the interface Σ.
For the present soluble case, however, the surfactant exists in the droplet domain Ωi as well as the interface Σ; thus we need to define two different concentrations (bulk concentration denoted by C(x, t) and interfacial concentration denoted by φ(x,t)) as described later. The fluid interface Σ separating the droplet and gas can be tracked in a Lagrangian manner X(α, t) = (R(α, t), Z(α, t)) with the parameter α ranging in [0, 1]. The unit tangent vector along the droplet interface is defined as τ = (τ1, τ2)= Xα/|Xα| with |Xα| = pR2α+ Z2α, where the subscript α indicates the partial derivative with respect to α. Therefore, the unit outward normal vector directing from the droplet Ωi into the gas Ωo can also be defined as n = (τ2,−τ1).
The immersed-boundary method (Peskin 2002) is a smoothing interface method that formulates the two-fluid problem into a unified one-fluid system by exploiting the indicator function Unverdi & Tryggvason (1992) (or Heaviside function):
H(x, t) = Z
Ωi
δ(x− x0)dx0. (3.1)
Here, H = 1 represents the droplet phase (with the viscosity µi and density ρi), while H = 0 represents the gas phase (with µo and ρo). Thus, this one-fluid system has
viscosity and density defined by µ = µiH + µo(1 − H) and ρ = ρiH + ρo(1 − H), respectively. Another advantage of using the indicator function is to simplify the numerical procedure for solving the bulk surfactant concentration. In the present flow with soluble surfactant, the bulk concentration C is defined only in the droplet phase Ωi but not in Ωo. One can redefine the bulk concentration by HC so that the convection–diffusion equation for the bulk concentration is solved in the regular computational domain Ω rather than in the irregular non-stationary droplet domain Ωi. Thus, we can regard HC as the effective bulk surfactant concentration. The surfactant absorption and desorption between the bulk and interface can be termed as a singular delta source along the interface in our immersed-boundary framework. The above formulation for two-phase flow with soluble surfactant in two dimensions (2-D) can be found in detail in the recent work of Chen & Lai (2014). In the next section, we describe the governing equations for interfacial flow with soluble surfactant in a 3-D axisymmetric formulation.
Before we proceed, we introduce the mathematical definitions for the gradient, divergence and Laplace operators in axisymmetric cylindrical coordinates:
∇= ∂
∂r, ∂
∂z
, ∇·=
1 r
∂
∂rr, ∂
∂z
, ∆= ∇ · ∇ =1 r
∂
∂r
r ∂
∂r
+ ∂2
∂z2. (3.2a−c) To compute the indicator function, one can first take the gradient operator for (3.1) and then apply the divergence operator to the resultant equation; so we obtain the following Poisson equation
1H(x, t) = −∇ ·Z
Σ
nδ(x− X) dΣ. (3.3)
The non-dimensional Navier–Stokes flow with soluble surfactant in the usual immersed-boundary formulation can be written as
ρ ∂u
∂t +(u · ∇)u
+ ∇p = 1 Re
∇ ·µ(∇u+ ∇uT)−µ r2u∗
+ f We +
ρg
Fr2, (3.4)
∇ ·u= 0, (3.5)
f= Z 1
0
F(α,t)δ(x − X(α, t)) dα, (3.6)
F(α,t) = ∂
∂α(σ (φ)τ(α, t)) − σ(φ)Zα
Rn(α,t), (3.7)
∂X(α,t)
∂t = U(α, t) = Z
Ω
u(x,t)δ(x − X(α, t)) dx, (3.8)
∂φ
∂t +(∇s·u)φ= 1 Pes
1 R|Xα|
∂
∂α
R
|Xα|
∂φ
∂α
+ Φ(α, t), (3.9)
∂(HC)
∂t + ∇ · (uHC) = 1
Pe∇ ·(H∇C) −Z 1 0
Φ(α,t)δ(x − X(α, t)) dα. (3.10) Under the assumption of axisymmetry, the velocity field u = (u, w) defined on the domain Ω has the radial u and axial component w. In (3.4), an additional term in terms of u∗ = (u, 0) is introduced into the radial momentum equation, which arises from the cylindrical coordinate system. Equation (3.7) describes the
immersed-boundary force arising from the surface tension σ which depends on the interfacial surfactant concentration φ. A well-known model, the Langmuir equation of state (Eastoe & Dalton 2000; Ceniceros 2003; Lin et al. 2003), has been used extensively as the constitutive law to specify the relationship between surface tension σ and surfactant concentration φ. However, the Langmuir model is not adequate for certain types of surfactant and can barely be used for low interfacial concentration.
In this work, we shall use the data fitting to the experimental measurements for the particular surfactant S111n to obtain the constitutive function σ (φ). The fluid interface, Σ, simply moves with the local fluid velocity as shown in (3.8). The interaction between the Eulerian and Lagrangian variables is rendered by the Dirac delta function δ(x) = δ(r)δ(z). The dimensionless numbers in the fluid equations are the Reynolds number (Re = ρiR0U∞/µi), the Weber number and the Froude number (Fr = U∞/√g∞R0). Here we define the characteristic scales by R0 the radius of droplet at rest, U∞ the relative impact speed, σ0 the clean surface tension and g∞ the gravitational constant.
Equations (3.9), (3.10) are the coupled convection–diffusion equations for the interfacial and bulk surfactant concentrations. In (3.9), ∂/∂t represents the material derivative since the interfacial surfactant φ(α, t) is defined in the Lagrangian coordinate α. The surface gradient operator is defined as ∇s = (I − nn)∇. Here, the coupled term Φ = SaCΣ(1 − φ) − Sdφ describes the absorption and desorption mechanisms between the interfacial and the bulk surfactant (Eggleton & Stebe 1998; Eastoe & Dalton 2000; Tabor, Eastoe & Dowding 2009). This coupling ensures that the total surfactant mass is conserved mathematically. Here, Sa and Sd are the adsorption and desorption constants, respectively, and CΣ is the bulk surfactant concentration evaluated adjacent to the interface. To set the bulk surfactant concentration as zero in the gas, we take advantage of the indicator function so that CΣ can be simply evaluated by
CΣ(α,t) = Z
Ω
HCδ(x − X(α, t)) dx. (3.11)
The dimensionless numbers for the surfactant equations are the bulk Peclet number Pe = U∞R0/Db and surface Peclet number Pes= U∞R0/Ds, where Db and Ds are the diffusivities of the bulk and interfacial surfactants, respectively. To make the present model well posed, we need proper boundary conditions for the velocity u and bulk concentration C on ∂Ω, and initial conditions for u(x, 0), C(x, 0), φ(α, 0) and the interface configuration X(α, 0).
The short-range forces are formulated based on the well-known DLVO theory (Israelachvili 2011), named after Derjaguin, Landau, Verwey and Overbeek, which accounts for the combined effects of van der Waals (vdW) attractive force and EDL repulsive force. The classical description of such intermolecular forces (Zhang et al.
2010; Conlisk et al. 2012) arises from the Poisson–Boltzmann equation. By treating the thin film between the surfaces of two colliding droplets as a layer separating two identical charged, parallel planes, specifically around the annular rim where there is the smallest thickness of air gap and the curvature is relatively small (Pan et al.
2008), one can express the intermolecular forces per unit volume as the gradient of the corresponding disjoining pressure and insert it simply as a source term in the IB formulation, which is a function of the distance between the two planes. An approximation of the disjoining pressure, pdj= pvdw+ pedl, is a combination of vdW effect, pvdw, and EDL effect, pedl, which are given in dimensional forms as
pvdw= − A0
6πh3g, pedl= E0exp(−κ0hg), E0= 64nkT tanh2z0eϕ 4kT
, (3.12a−c)
where hg(r) is the width of the gap at r, A0 is the Hamaker constant (typically 10−21 to 10−19 J for water), E0 is the Debye constant, κ0−1 is the Debye length (960 nm for water in liquid, which is approximately eight times smaller in air), k is the Boltzmann constant, T is the absolute temperature, z0 is the charge number, e is the elementary electron charge, n is the ion concentration and ϕ is the surface potential (Zhang et al. 2010). These two intermolecular factors affect the fluid-flow system through a volumetric force field in the dimensionless form as
fmolecule=∂pdj
∂r = A
ˆh4g − Eκ exp(−κ ˆhg)
!∂ ˆhg
∂r . (3.13)
Here the dimensionless Hamaker constant is A = A0/(2πρiR30U2∞) and the Debye constant is E = E0/(ρiU2∞). κ and ˆhg are normalized variables.
3.2. Numerical algorithm
In the numerical processes, all the fluid variables are defined at the cell centre labelled as xij= (ri−1/2,zj−1/2)= ((i − 1/2)h, c + ( j − 1/2)h) in Ω with the uniform mesh width h in the r and z directions. It is noted that the uniform spatial grid used here is just for presentation purposes and it can be extended to a stretching grid without great difficulty in practice. For the immersed interface, we use a collection of discrete Lagrangian markers Xk = X(αk)= (Rk,Zk) to track the interface where αk= k1α, k = 0, 1, . . . , M are the parametric points. Typically, the size of 1α is chosen as 1α ≈ h/2. The discrete value fk denotes an approximation evaluated at αk, while fk+1/2 denotes the approximation evaluated at αk+1/2= (k + 1/2)1α. Using the standard centred-difference approximation, the unit tangent τk+1/2, the unit normal nk+1/2 and the stretching factor |Xα|k+1/2 are all evaluated at the half-integer index points. In addition, we also define the interfacial surfactant concentration φk+1/2 and the surface tension σk+1/2 at those points with half-integer indices.
Let 1t be the size of the time step, and the superscript n be the index. At the beginning of each time step, e.g. step n, the interface configuration Xnk, the fluid velocity unij, the interfacial surfactant concentration φnk+1/2 and the bulk concentration Cnij must be given. Hence the complete numerical algorithm for solving (3.4)–(3.10) to advance one time step can be summarized as follows.
(1) Compute the surface tension σk+1/2n and unit tangent τk+1/2n , and then calculate the interface force Fnk as in (3.7).
(2) Distribute the interfacial force from the Lagrangian markers into the fluid as in (3.6) to obtain fnij.
(3) Solve the Navier–Stokes equations (3.4), (3.5) by the projection method to update the new fluid velocity un+1ij .
(4) Interpolate the new velocity on the fluid lattice points onto the marker point and move the marker point to update interface position Xn+1k as shown in (3.8).
(5) Solve the equation (3.3) to obtain the new indicator function Hijn+1.
(6) Solve the equation of interfacial surfactant concentration (3.9) to update φn+1k+1/2. (7) Solve the bulk concentration equation (3.10) to update Cijn+1.
The above steps (1–4) are the conventional numerical procedures for the immersed- boundary method, for which the details are not reiterated. Readers who are interested in the numerical implementation in 2-D and 3-D axisymmetric Navier–Stokes flows
with insoluble surfactant can find details in Lai et al. (2008), Lai, Huang & Huang (2011). Here, we focus on the numerical steps (5–6) since the treatment of soluble surfactant is the main numerical topic of this paper.
In step (5), the indicator function can be discretized via the Poisson equation (3.3) directly using the centred-difference scheme by
∆hHijn+1= −∇h· X
k
nn+1k δh(xij− Xn+1k )|Xα|n+1k 1α
!
, (3.14)
where ∆h and ∇h· are the standard centred-difference approximations for the Laplace and the divergence operators in axisymmetric cylindrical coordinates, respectively. The discrete delta function δh is adopted as the smoother one to avoid the oscillation of immersed force calculations, as introduced in Yang et al. (2009). The normal vector and the stretching factor evaluated at the marker Xk can be simply averaged by nk= (nk+1/2 + nk−1/2)/2 and |Xα|k = (|Xα|k+1/2 + |Xα|k−1/2)/2, respectively. The accuracy of the above numerical computation for the indicator function is well investigated in Chen et al. (2011).
To update the interfacial surfactant in step (6), we first rewrite the interfacial surfactant equation by multiplying the surface stretching factor R|Xα| on both sides of (3.9). We then substitute the surface differential relation ∂/∂t(R|Xα|) = (∇s·u)R|Xα| Lai et al. (2011) into the resultant equation to obtain
∂φ
∂tR|Xα| +∂(R|Xα|)
∂t φ= 1 Pes
∂
∂α
R
|Xα|
∂φ
∂α
+ Φ(α, t)R|Xα|. (3.15) By applying the implicit Euler method for the time integration and using the centred- difference scheme for the spatial discretization, we have
(φR|Xα|)n+1k+1/2− (φR|Xα|)nk+1/2
1t − Φk+1/2n+1 (R|Xα|)n+1k+1/2
= 1
Pes1α
Rn+1k+1
|Xα|n+1k+1
φk+3/2n+1 − φk+1/2n+1
1α − Rn+1k
|Xα|n+1k
φk+1/2n+1 − φk−1/2n+1
1α
!
, (3.16)
Φk+1/2n+1 = SaC∗k+1/2(1 − φn+1k+1/2)− Sdφk+1/2n+1 , (3.17) C∗k+1/2=X
ij
Hijn+1Cnijδh(xij− Xn+1k+1/2)h2. (3.18)
Although the above scheme seems complicated at first glance, this discretization results in a tridiagonal symmetric linear system for φk+1/2n+1 which can be solved very easily.
As for solving the bulk surfactant in step (7), we discretize (3.10) in a similar manner (implicit Euler method in time and centred-difference method in space) as in step (6) to obtain
(HC)n+1ij − (HC)nij
1t + ∇h·(uHC)nij= 1
Pe∇h·(H∇hC)n+1ij − ˜Cijn+1, (3.19)
˜Cn+1ij =X
k
Φk+1/2n+1 δh(xij− Xn+1k+1/2)1α. (3.20)
Since the coupling terms Φk+1/2n+1 for all k are already known from step (6), the above discretization involves solving a variable diffusion equation for the bulk concentration Cn+1ij . It is noted that, in practical numerical implementation, we need to regularize the diffusion coefficient (indicator function H) by using √
H2+ 2 instead of H simply because H equals to zero in the gas phase Ωo. Here, is chosen to be approximately 10−6 as suggested in Chen & Lai (2014).
3.3. Validation
The simulation is based on the similar methodology, i.e. the immersed-boundary method and front tracking, which has been applied in our previous studies for droplet related problems (Pan & Law2007; Pan et al.2008). Substantial agreement regarding the evolution of deformed surfaces between the experiment and numerical simulation has been demonstrated, as also shown in the appendix A. Moreover, the advanced algorithm including surfactant dynamics is included in Chen & Lai (2014), in which a reasonable procedure of validation is performed. In addition, a convergence test along with the conservation of soluble surfactant is to be provided in §5.2 with more discussion.
4. Experimental observation on impact outcome with addition of surfactant The typical scenarios of droplet collision are shown in figure 2. Here the droplets of water generally merge at low/moderate Weber numbers and impact parameters (figure 2a,b). Bouncing can only occur when the colliding path is significantly off centred, as shown in figure 2(c). At higher We, after merging, separation of the coalesced droplets is created in order to balance the excess energy that cannot be contained in a single droplet. Stretching and reflexive (near head-on) forms of separation are observed in large (figure 2d) and small B (figure2e), respectively. With addition of surfactant, however, bouncing is generated even in head-on collision.
Figure 5 shows the regime diagram of droplet impacts. It is noted that, in the experiments, dimensional values such as surface tension are recorded. They will however be transformed to dimensionless quantities in order to be compared with computational results as discussed later. For pure water (figure 5a), merging is the typical outcome right after collision, whether it is permanent at low We or temporary coalescence at high We, followed by separation of the merged droplet and further generation of satellite droplets. Bouncing of the impinging drops occurs only when B is so large that the gas film remains in between and may act as a buffer to rebound the droplets. This regime was generally overlooked in early studies (Brazier-Smith et al. 1972; Ashgriz & Poo 1990; Jiang et al. 1992), but was delineated in the study of Qian & Law (1997). While nitrogen was used as the environmental gas, the regime reported therein was similar to the present one for which air was adopted, whose properties deviates little from that of nitrogen.
By adding a small amount of surfactant, approximately 0.1 % (the initial mass fraction in a droplet, designated by ψ0), however, the bouncing regime (II) is substantially enlarged, as shown in figure 5(b). In particular, even in head-on collisions with B = 0, the droplets rebound after impact, which has never been observed for water. When more surfactant is added, as shown in figure 5(c) for a 0.3 % concentration of S111n, which is slightly below the CMC limit (see figure4a), the bouncing regime shows significant expansion. For a concentration far beyond the CMC limit, nevertheless, this regime is ostensibly shrunk, as shown in figure 5(d) for 1 % S111n.
0 5 10 15 20 25 30 35 40 0.1
0.2 0.3 0.4 B 0.5
B
We We
0.6 0.7 0.8 0.9
(a) 1.0 (b)
(c) (d)
(II)
(II)
(II)
(II) (V)
(V) (V)
(III)
(III)
(III) (I)
(I)
(V)
(IV) (IV)
(IV) (IV)
0 10 20 30 40 50 60 70 80 90 100 110
0 10 20 30 40 50 60 70 80 90 100 110 0 10 20 30 40 50 60 70 80 90 100 110
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
FIGURE 5. (Colour online) The regime diagram of impact on (We, B) for (a) pure water with σ0= 0.072 N m−1; (b) ψ0= 0.1 % surfactant S111n with σ0= 0.032 N m−1; (c) ψ0 = 0.3 % with σ0= 0.021 N m−1; (d) ψ0= 1.0 % with σ0 = 0.017 N m−1. The diameter of droplets was fixed at ∼300 µm except for pure water, which ranged from 280 to 440 µm. The impact velocity varied from 0.2 to 2.6 m s−1. The symbols in the diagrams refer to the regimes as indicated in figure 1; here open circles indicate coalescence, stars bouncing and crosses separation. In (a), since the coalescence regimes at low B are not intervened by the bouncing regime (II) created only at high B, they are designated together by I/III.
Although the reduction in surface tension may play a role in such variations, as observed for other natural liquids such as hydrocarbons (Jiang et al. 1992), it is not the major cause for the present cases with a surfactant. For demonstration, we have also adopted different surfactants (S131 and S386) while keeping surface tension fixed, so the mixture of water includes much less surfactant, i.e. 0.005 %, as shown in figure 6(a) for S386 that is non-ionic. In contrast to figure 5(b) for S111n, the bouncing regime is not enlarged much with the addition of these surfactants. However, when more surfactant is added (figure 6b), specifically when approaching the CMC limit, remarkable expansion of the bouncing regime is again observed.
The experimental results reveal a predominant influence of the surfactant concent- ration, in addition to surface tension and the type of surfactant, which is related to the molecular structure of the chemical formulation and the effectiveness in reduction of surface energy. Such enhancement of rebounding between two aqueous droplets with surfactant can be understood by using numerical simulation for the correlations between the dynamics of the gas gap and the motions of fluids and surfactant.
5. Computational analyses and discussion of observed results
Figure 7(a) illustrates the geometry of the collision between two droplets (one from the top and the other from the bottom), where dc and dr are indicated in a
B
We We
(a) (b)
(II)
(II) (V)
(V)
(III) (IV) (IV)
0 10 20 30 40 50 60 70 80 90 100 110 0 10 20 30 40 50 60 70 80 90 100110
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
FIGURE 6. (Colour online) The regime diagram of impact on (We, B) for (a) ψ0= 0.005 % surfactant S386 with σ0 = 0.032 N m−1; and (b) ψ0 = 1.0 % with σ0 = 0.021 N m−1. The diameter of droplets was fixed at ∼300 µm. The impact velocity varied from 0.2 to 2.6 m s−1. The symbols in the diagrams refer to the regimes as indicated in figure 1.
z
0 0.002 0.004 0.006 0.008 0.010
(a) (b) (c)
r
1.2 0.8 0.4 0 -0.4 –0.8 –1.2
4h 2h h
t
0.5 1.0 1.5 2.0 2.5 3.0 0
0.005 0.010 0.015
Insoluble Soluble
Insoluble Soluble
Insoluble Soluble
t
1 2 3
Mass
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Arclength
0.5 1.0 1.5 2.0 2.5 3.0 0
0.1 0.2 0.3 0.4 0.5
Arclength
0 0.5 1.0 1.5 2.0 –0.4
–0.3 –0.2 –0.10 0.1 0.2 0.3 0.4 Center
Rim
(d) (e) ( f )
FIGURE 7. (Colour online) Convergence test and the effect of surfactant solubility: Re = 213.6, A = 1.4255 × 10−11, We = 8.9566 and Pes = 100. The simulations with soluble surfactant consider additional parameters, Pe = 100 Sa= 2.32 and Sd= 20. (a) Illustration of the deformed surfaces during impact, showing an exaggerated view of the gap between the droplets. (b) Grid refinement test based on the geometry and width of the gap at t = 1.0. (c) Evolution of dr for insoluble surfactant and soluble surfactant. (d) Evolution of surfactant mass. Here Mb indicates the mass in the bulk fluid, Ms is that on the surface and Mt= Mb+ Ms is the total mass in the system. (e) Distribution of interfacial surfactant concentration, φ, for insoluble (dash line) and soluble (solid line) cases at t = 0.5. ( f ) Distribution of Marangoni stress for insoluble (dash line) and soluble (solid line) cases at t = 0.5.
blow-up of the impact region. They are used to show the geometrical evolution of the gap expressed by the separation distance between the impinging surfaces (near the central plane which has a symmetric condition in the computation domain) of the two droplets at the centre and at the rim. It is noted that dr is measured at the radius where the minimum gap occurs, which is time dependent and is not
measured at a fixed radius. As interpreted in Pan et al. (2008), the surfaces become indented when the droplets approach each other and pressure increases in the gap. The impedance due to the gap pressure is critical for the occurrence of bouncing, which is strongly related to the drainage of the trapped gas, and merging can happen only if the distance between the interfaces is reduced to such an extent that short-range factors are effective. Therefore, the smallest width of the gap computed based on the Navier–Stokes continuum mechanics could be identified as an index of tendency toward bouncing. In particular, if the minimum gap width that can be attained during the collision period, typically at the rim (and hence the smallest dr, as designated by dr,min), is much larger than the scale below which intermolecular attractive forces or other short-range mechanisms become dominant, the repulsive pressure yields rebounding of the droplets. The minimum gap width is changed by the addition of surfactant, which may induce the effects of reduced surface tension, non-uniform distributions of surfactant and repulsive EDL forces. They are investigated herein computationally based on the numerical simulation.
In this section, numerical parameters are selected from the experimental conditions.
Instead of using dimensional quantities recorded for fluid properties in the experiments, the variables are non-dimensionalized in the presentation of computational results.
First of all, a convergence test and a study of the capability of the numerical method in capturing the essential dynamics of surfactant-laden fluids are performed. Secondly, the concentration of surfactant is varied to identify the trend and the cause that prevents droplets from merging in collisions as observed in the experiments. A series of numerical tests related to the experimental situations of head-on collisions between droplets are conducted, including an investigation of the difference in using insoluble and soluble surfactants, which has been rare in the literature, and clarification of the surface tension effect. By means of such an a posteriori discussion, it is shown clearly that the key mechanism can be related to the Marangoni stress. Furthermore, consideration of a wide range in surfactant amount shows that the enhancement of the bouncing tendency has a limitation owing to the CMC limit of surfactant concentration. In addition, to understand the role of the other properties of the liquid played in such variations, specifically of viscosity, water droplets are mixed with glycerol to contrast the factor and show whether the observed bouncing regime due to the presence of surfactant is affected by viscosity.
5.1. Numerical set-up
The conditions of numerical tests are based on those of the experiments, where the radius of a spherical water droplet is R0 = 147.5 µm and the impact speed is accelerated from zero by the gravitational field up to U∞ = 1.48 m s−1. The density and viscosity of the droplets, respectively, are ρi = 998 kg m−3 and µi = 1020 µPa s. The surface tension σ0 of pure water in 1-atm regular air is approximately 0.072 N m−1, and the density and viscosity of air are respectively, ρo = 1.2 kg m−3 and µo = 18.6 µPa s. Accordingly, the surface area, volume and mass of a droplet are Ad = 2.7340 × 10−7 m2, Vd = 1.3442 × 10−11 m3 and Md= 1.3415 × 10−8 kg, respectively.
These physical quantities are selected as the characteristic scales for normalization of the governing equations, and the relevant dimensionless numbers are Re = 213.6, We = 8.9566 and A = 1.4255 × 10−11. As discussed in Stebe, Lin & Maldarelli (1991), since the usual circumstance of mass transfer in liquids is characterized by high Peclet numbers, here we choose Pe = Pes= 100 just for simplicity, as in Eggleton & Stebe
0 0.5 1.0 1.5 2.0 2.5 0.1
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Experimental data
Monotonic exponential model
FIGURE 8. (Colour online) Curve fitting for the experimental data of S111n in figure 4.
The solid line is formulated by a monotonic exponential model and the diamond markers indicate the non-dimensionalized data of S111n.
(1998) and Valkovska & Danov (2000). In addition, the selection of Sa= 2.32 and Sd= 20 are within the range of consideration in Zhang, Eckmann & Ayyaswamy (2006).
Throughout this section, a computational domain Ω = [0, 3] × [0, 3] is used in all numerical tests.
The experimental data of surface tension versus surfactant concentration were illustrated earlier in figure 4. Table 2 lists the variation of surface tension with surfactant (S111n) concentration in a static droplet. The amount of surfactant ψ0 is measured as the percentage of its mass, Md, initially contained in a droplet.
Figure 8 depicts the relation of the dimensionless surface tension σ and interfacial concentration φ (for S111n) in terms of a polynomial fitting curve. It is seen that surface tension decreases dramatically in the beginning and then tends to level off after ψ∞= 0.35 % (see figure 4) which we name it as the maximum surfactant package.
In this way, the dimensionless interfacial concentration of surfactant φ = ψ0/ψ∞ is used for a given ψ0 in the numerical processes and it is normalized to unity when CMC is attained. For instance, the mass concentration ψ0= 0.1 % gives the value of φ= 0.286.
According to table 2, the constitutive equation of σ -φ can be approximated by a monotonic exponential model (dimensionless formulation):
σ (φ)=
0.754e−3.45φ+ 0.246, φ < 1
0.270, φ> 1. (5.1)
Since the data in table 2 were measured at steady state, we simply assume that all the surfactant already attached to the interface when ψ0< ψ∞, and the corresponding φ is simply ψ0/ψ∞. If ψ0> ψ∞, the surface has reached the maximum package filled with surfactant and the remainder are dissolved into the droplet. For instance, if we add 0.5 % surfactant into the droplet, then φ0.5 %= 1 and C0.5 %= 3/7.
5.2. Convergence test, solubility and mass conservation
To understand evolution of the intervening gap during droplet collision, we define two characteristic distances; namely, dc (gap at the centre) and dr (gap at the rim), which are used to present the gap dynamics. Figure 7(a) shows the impact diagram of two droplets. The code is validated through a mesh refinement process, and the
ψ0(%)00.10.110.20.220.30.390.40.440.50.560.710.80.9 σ0(mNm−1)723230.82524211720.5171916.617.21717.1 |eσ|00.1860.1640.0310.0090.040.0010.1420.0530.1420.0220.1290.1420.136 TABLE2.Experimentalmeasurementofsurfacetensionσ0(mNm−1 )versusS111nsurfactantfraction(ψ0,masspercentageofMd).The relativeerrorofthemonotonicexponentialfittingaccordingto(5.1),|eσ|,isalsoshown.Here|eσ|isdefinedas(σexp−σfit)/σexp,whereσexp isthedimensionlesssurfacetensionobtainedfromtheexperimentwhileσfitisthatfromthemonotonicexponentialmodel.σfitisfixedat 0.754e−3.45 +0.246asψ0>0.35%.
convergence of the numerical method is explained by the profile of droplet surface at t = 1.0 which is approximately the time when the minimum gap is attained.
Figure 7(b) shows three profiles of the deforming gap with minimum mesh sizes 4h, 2h and h, respectively, where a linear convergence is observed. It is noted that, with variation of the gap size, dr, at least two cells are always kept in between two interfaces, particularly at the rim, such that the IB calculation is valid.
To gain an insight into the effect of surfactant solubility on the dynamics of droplet collision we consider two cases, both having ψ0 = 0.1 % (so φ = 0.286) surfactant inserted but one is soluble and the other is insoluble. Figure 7(c) shows the gap dynamics before t = 3.0. The dashed line represents the evolution of dr in the insoluble case and the solid line represents the result of the soluble case. Obviously, the insoluble case has a bigger gap size at the rim, as compared to that of the soluble case. A reasonable explanation is that the insoluble case keeps all surfactant on the interface such that a stronger Marangoni stress (compared to the soluble case) due to non-uniform distribution of surfactant concentration suppresses more the air drainage during the collision (to be interpreted further in the next section). The corresponding distributions of surfactant concentration and Marangoni stress at t = 0.5 are shown in figure 7(e,f ), respectively. It is seen that large amplitudes appear in both cases outside the rims of the deformed droplet surfaces, which are indicated by circular markers.
Specifically, the negative stresses closer to the rims indicate the inward direction of surfactant induced forces that would prevent the air flow in the gap from draining out, as will be further demonstrated.
To further test the numerical accuracy and solubility of surfactant, we have examined the natural oscillation of a droplet for cases with insoluble and soluble surfactant. The oscillation period is predicted for a clean droplet as T0= 2π(ρiR30/8σ0)1/2 (Lamb 1932), based on an inviscid condition with small oscillation amplitude. With surface tension set identically to 0.032 N m−1 at the beginning, when surfactant is uniformly distributed in the droplet for both cases, the oscillation periods are the same, i.e. T0∼ 7.0277 × 10−4 s. As time passes however, the soluble case exhibits a shorter duration, T0= 5.9797 × 10−4 s. This can be reasonably explained by the reduction of interfacial concentration, hence increase of surface tension, due to the solubility of surfactant.
Figure 7(d) shows the evolutions of Mb (the surfactant mass in the bulk fluid), Ms
(the surfactant mass on the interface) and the total mass Mt of the system. Due to the solubility of surfactant, the mass on the interface Ms (dashed line) decreases with time while the mass in the bulk fluid Mb (dotted line) increases. Since the mass-preserving numerical scheme guarantees perfect conservation of surfactant mass with minimum roundoff error in a discrete sense through summation by parts, the total mass Mt
remains invariant for both soluble and insoluble cases, and the two curves of evolution overlap.
5.3. Marangoni effect due to surfactant
In order to clarify the key factors that cause the water droplets to bounce off instead of coalescence in the presence of surfactant, we first focus on the magnitude of dimensionless surface tension, σ . Variation of surface tension along the droplet interfaces may induce effects of both capillary stress and Marangoni stress. To identify these using computational analysis, we test two cases for clean droplets (no surfactant), i.e. one with σ = 1 and the other with σ = 0.41. The effects can be demonstrated via the variations in the neighbourhood of the thin air film between the
