1. Introduction
Small living creatures such as water striders, beetles and mosquito larvae use surface tension to stand, walk, leap or hang on the surface of water (Hu & Bush Reference Hu and Bush2005; Bush & Hu Reference Bush and Hu2006; Feng et al. Reference Feng, Gao, Wu, Jiang and Zheng2007; Vella Reference Vella2015; Lee, Kim & Lee Reference Lee, Kim and Lee2017). Inspired by these natural occurrences, researchers have developed millimetre-scale robots (Koh et al. Reference Koh, Yang, Jung, Jung, Son, Lee, Jablonski, Wood, Kim and Cho2015; Hu et al. Reference Hu, Lum, Mastrangeli and Sitti2018) for transport across the surface of a liquid that might be useful in targeted drug delivery, minimal invasive surgery, and other bio-engineering applications. These robots feature a hydrophobic surface with strong interfacial tension that prevents the body from breaking the liquid surface and sinking. Once the body rests at the surface, additional locomotion can be provided utilizing the techniques described by Hu et al. (Reference Hu, Lum, Mastrangeli and Sitti2018), Jiang et al. (Reference Jiang, Gao, Zhang, He, Zhang, Daniel and Yao2019) and Grosjean, Hubert & Vandewalle (Reference Grosjean, Hubert and Vandewalle2018). Many other biomedical applications require encapsulation of one liquid in another. Examples include separation (Albertsson Reference Albertsson1986; Zhang et al. Reference Zhang, Cai, Lienemann, Rossow, Polenz, Vallmajo-Martin, Ehrbar, Na, Mooney and Weitz2016; Li et al. Reference Li, Xie, Liu, Kong, Song, Wen and Jiang2018) or encapsulation (Orive et al. Reference Orive2003; Delcea, Möhwald & Skirtach Reference Delcea, Möhwald and Skirtach2011) of biomolecules and cells. In this context, aqueous two-phase systems (Hann et al. Reference Hann, Niepa, Stebe and Lee2016; Hann, Stebe & Lee Reference Hann, Stebe and Lee2017; Chao et al. Reference Chao, Mak, Rahman, Zhu and Shum2018; Xie et al. Reference Xie, Forth, Chai, Ashby, Helms and Russell2019), formed using a mixture of dextran and poly(ethylene glycol) (PEG) which phase separates to form two immiscible aqueous phases, are widely used.
We perform three-dimensional numerical simulations on two immiscible aqueous fluids to demonstrate that a droplet of higher density can either hang from the surface like mosquito larvae, bounce on the surface like water striders, or form a shroud that completely wraps the denser fluid as it sinks in the pool of a lighter liquid. As the drop makes contact with the pool, the evolution of the three-phase contact line (TPCL) plays a major role in the dynamics of a drop hanging or sinking from the surface. It will be shown using force balance equations that during the hanging process, the surface tension force balances the heavier droplet at the surface of the pool. The size of the droplet and its initial kinetic energy are two of the key parameters that dictate the outcome in this situation. Xie et al. (Reference Xie, Forth, Zhu, Helms, Ashby, Shum and Russell2020) presented experimentally a similar phenomenon of hanging (Phan et al. Reference Phan, Allen, Peters, Le and Tade2012; Phan Reference Phan2014) and wrapping (Kumar et al. Reference Kumar, Paulsen, Russell and Menon2018) using an aqueous two-phase system of a dextran solution containing polycations, and a PEG solution containing polyanions. In the presence of oppositely charged polyelectrolytes, after coming into contact, the solutions create structured coacervate sacs of negligible mass and thickness at their interface. These coacervate sacs effectively increase the interfacial tension between the two solutions, resulting in hanging of the heavier droplets from the pool surface.
2. Methods
The volume of fluid (VOF) approach of Hirt & Nichols (Reference Hirt and Nichols1981) serves as a foundation for calculations involving two fluids separated by a sharp interface. The VOF approach achieves excellent compliance with mass conservation, but it can be difficult to capture the geometric features of a complex interface. Osher & Sethian (Reference Osher and Sethian1988) introduced the level set (LS) method, which is an efficient interface capture technique. This approach captures the interface properly, although it may violate mass conservation in some circumstances. A combination of the LS approach with the VOF method, known as the coupled level set and volume of fluid (CLSVOF) method, can accomplish mass conservation and capture the interface properly. The LS function is utilized exclusively to compute the geometric characteristics at the interface in the CLSVOF technique (Sussman & Puckett Reference Sussman and Puckett2000), while the volume fraction is determined using the VOF method. Continuum surface tension force by Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992) has been used widely to evaluate the source term due to surface tension. However, a free-energy-based surface tension force model is proposed by Yuan et al. (Reference Yuan, Chen, Shu, Wang, Niu and Shu2017) for simulation of multi-phase flows by the LS method, which outperforms the previous continuum surface tension force model in terms of accuracy, stability, convergence speed and mass conservation. Howard & Tartakovsky (Reference Howard and Tartakovsky2021) also extended the conservative LS method for 
$N$ fluid phases by introducing a new compression–diffusion equation that handles large deformation and triple junctions more accurately. In order to solve the 
$N$-phase flow problems, algorithms with 
$N$ (Ruuth Reference Ruuth1998), 
$N-1$ (Smith, Solis & Chopp Reference Smith, Solis and Chopp2002; Zlotnik & Díez Reference Zlotnik and Díez2009), 
$N(N-1)/2$ (Starinshak, Karni & Roe Reference Starinshak, Karni and Roe2014a,Reference Starinshak, Karni and Roeb) and 
$\log _2N$ (Chan & Vese Reference Chan and Vese2001) LS functions have been used.
2.1. Governing equations
Considering incompressible Newtonian fluids, the mass and momentum conservation equations for fluids 1, 2 and 3 are given by
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U} = 0, \end{gather}
\begin{gather}\rho\left(\frac{\partial \boldsymbol{U}}{\partial t} + \boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{U}\right) ={-}\boldsymbol{\nabla}\mathcal{P} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(2\mu\boldsymbol{D} \right) + \boldsymbol{F} + \boldsymbol{F_{st}}, \end{gather}
where 
$\boldsymbol {U}$ is the velocity vector field with components 
$({U}_1, {U}_2, {U}_3)$, 
$\mathcal {P}$ represents the dynamic pressure, 
$\rho$ and 
$\mu$ are scalar fields representing density and dynamic viscosity, 
$\boldsymbol {D}$ is the deformation tensor, and 
$\boldsymbol {F}$ and 
$\boldsymbol {F_{st}}$ are body force and surface tension force per unit volume. Here,
\begin{equation} \boldsymbol{D} = \tfrac{1}{2}\left( \boldsymbol{\nabla} U + \boldsymbol{\nabla} U^{\rm T}\right). \end{equation}Gravitational force is the only body force acting on all the fluids in our case. Surface tension force, as given by Howard & Tartakovsky (Reference Howard and Tartakovsky2021), is used in the momentum equation as
\begin{gather} \rho\left(\frac{\partial \boldsymbol{U}}{\partial t} + \boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{U}\right) ={-}\boldsymbol{\nabla}\mathcal{P} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(2\mu\boldsymbol{D} \right) + \rho\boldsymbol{g} + \sum_{i,j=1}^{3} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{3}{2}\,\sigma_{ij}\varepsilon\, \boldsymbol{\nabla}\varphi_i\times\boldsymbol{\nabla}\varphi_j\right). \end{gather}
Here, 
$\boldsymbol {g}$ is the acceleration due to gravity, and 
$\sigma _{ij}$ represents the surface tension at the interface between the fluids 
$i$ and 
$j$. The interfacial numerical thickness 
$\varepsilon$ is defined based on grid size 
$\Delta x$ as 
$\varepsilon =k_\varepsilon \,\Delta x$; for all the simulations reported here, 
$k_\varepsilon =1.5$ is used. The surface tension force is obtained using the free energy surface tension force model. The free energy density for 
$N$ immiscible fluids is given by Dong (Reference Dong2014).
In the present work, CLSVOF is used, which combines the advantages of both the LS method and the VOF method. The LS function is defined as a signed distance function from the phase interface such that
\begin{equation} \phi_i(\boldsymbol{x}) \left\{ \begin{array}{@{}ll} > 0 & \mbox{if }\boldsymbol{x}\in\varOmega_i, \\ = 0 & \mbox{if }\boldsymbol{x}\in\varGamma_i, \\ < 0 & \mbox{if }\boldsymbol{x}\notin\varOmega_i, \end{array} \right. \end{equation}
where 
$\varOmega _i$ denotes the subdomain containing the fluid of the 
$i$th phase, and 
$\varGamma _i$ is the sharp interface of the 
$i$th phase. The VOF function 
$f_i$ is the fraction of the cell volume occupied by the fluid of phase 
$i$ and satisfies the condition
\begin{equation} \sum_{i=1}^{3}f_i = 1. \end{equation}
The scalar field 
$\varphi _i$ used in (2.4) is defined using the LS function as
\begin{equation} \varphi_i = {H}(\phi_i). \end{equation}Here, the Heaviside function is defined as
\begin{equation} {H}(\phi) = \left\{ \begin{array}{@{}ll} 0 & \mbox{if} \ \phi <{-} \varepsilon, \\[4pt] \dfrac{1}{2}\left [ 1+\dfrac{\phi }{\varepsilon} + \dfrac{1}{{\rm \pi} }\sin \left (\dfrac{{\rm \pi} \phi}{\varepsilon } \right ) \right ] & \mbox{if} \ \left | \phi \right | = \varepsilon, \\[9pt] 1 & \mbox{if} \ \phi >{+} \varepsilon. \end{array} \right. \end{equation}The varying density and viscosity fields are also defined using the Heaviside function as
\begin{gather} \rho = \rho_1\,{H}(\phi_1)+\rho_2\,{H}(\phi_2)+\rho_3\,[1-{H}(\phi_1)-{H}(\phi_2)], \end{gather}
\begin{gather}\mu = \mu_1\,{H}(\phi_1)+\mu_2\,{H}(\phi_2)+\mu_3\,[1-{H}(\phi_1)-{H}(\phi_2)]. \end{gather}The motion of interfaces is tracked by solving explicitly the advection equation for both LS and VOF functions:
\begin{gather} \frac{\partial \phi_i}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{U}\phi_i) = 0, \end{gather}
\begin{gather}\frac{\partial f_i}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{U} f_i) = 0. \end{gather}2.2. Boundary conditions
The governing equations are solved in a three-dimensional Cartesian space. A closed system is considered for the simulations such that no fluid enters or leaves the computational domain:
\begin{equation} \boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{n} = 0. \end{equation}The no-slip boundary condition is assumed at all the boundaries of the computational domain. The boundaries of the computational domain are kept sufficiently away from the droplet to ensure that it does not affect the dynamics of the flow:
\begin{equation} \boldsymbol{n}\times\boldsymbol{U} = 0. \end{equation}Therefore, a Dirichlet boundary condition is used for the velocity field on all the boundaries. On the other hand, a Neumann boundary condition is used for pressure at the boundaries:
\begin{equation} \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla}\mathcal{P} = 0. \end{equation}2.3. Numerical methods
The governing partial differential equations are advanced in time using an explicit third-order Runge–Kutta method (Williamson Reference Williamson1980). A staggered grid is used for the discretization in space where the vector field quantities (such as 
$\boldsymbol {U}$) are defined at the cell face centre and the scalar field quantities (
$\mathcal {P},\rho,\mu,\phi,f$) are defined at the cell centre. The advection terms are discretized using a second-order ENO scheme as used by Chang et al. (Reference Chang, Hou, Merriman and Osher1996) and Son & Dhir (Reference Son and Dhir2007). The viscous terms are discretized using a second-order central difference scheme. We note here that the viscosity 
$\mu$ is not constant throughout the domain, so special care has to be taken to include 
$\mu$ in the discretization scheme. The pressure Poisson equation, which is employed to project a velocity field into a divergence-free space, is solved using a parallel multigrid iterative solver (Pal Reference Pal2020; Pal & Chalamalla Reference Pal and Chalamalla2020) to obtain the dynamic pressure. To advance in time for the advection equations of the LS and VOF functions, we use an operator splitting algorithm (Son Reference Son2003), in which we solve (2.11) and (2.12) in one direction at a time. The operator splitting is of second-order accuracy in time, and the order of sweep direction at each time step is also alternated. The solution of the advection equation for the LS function does not satisfy the signed distance property from the interface. For this, the LS function is re-initialized at each time step after the operator splitting algorithm (Son Reference Son2003).
To perform a three-phase flow simulation (
${N}=3$), only two (
${N}-1$) phase equations are solved using the CLSVOF algorithm. The numerical solution to these 
${N}-1$ phase equations generates some voids and overlaps between the phases. By using the 
${N}-1$ phase equation, it is assumed that the 
${N}$th phase occupies the void region, and this avoids the singularity problems in the computational domain:
\begin{gather} f_3 = 1-f_1-f_2, \end{gather}
\begin{gather}{H}(\phi_3)=1-{H}(\phi_1)-{H}(\phi_2). \end{gather}A VOF correction is performed to overcome the overlap issues, such that
\begin{equation} f_2 = 1 - f_1 \quad \mbox{if } f_1+f_2> 1. \end{equation}This VOF correction is biased towards phase 1 as we assume that phase 1 represents the primary fluid of interest. The coupled nature of the CLSVOF algorithm appropriately adjusts the LS function for this VOF correction.
A constant time step size is used such that it satisfies the following time step restrictions. First, the standard Courant–Friedrichs–Lewy (CFL) condition is satisfied:
\begin{equation} \Delta t_u \leqslant \frac{CFL}{\dfrac{|u|_{{max}}}{\Delta x} + \dfrac{|v|_{{max}}}{\Delta y} + \dfrac{|w|_{{max}}}{\Delta z}}. \end{equation}According to Brackbill et al. (Reference Brackbill, Kothe and Zemach1992), when treating the surface tension term explicitly, the time step must be sufficiently small to resolve the capillary waves phenomena. This gives another time step restriction as
\begin{equation} \Delta t_\sigma \leqslant CFL_\sigma \sqrt{\frac{\min(\rho_i+\rho_j)\times\min(\Delta x,\Delta y, \Delta z)^3}{\max(4{\rm \pi}\sigma_{ij})}}, \quad i\ne j. \end{equation} We have used 
$CFL=0.5$ and 
$CFL_\sigma = 0.5$ for all cases. Other time step restriction criteria based on viscosity and gravity give more relaxed values. A constant time step is used such that it satisfies the above-mentioned restrictions sufficiently throughout the simulation.
3. Validation of numerical approach
3.1. Advection test
In this study, a parallel three-phase incompressible flow solver is used that is an extension of an existing two-phase flow solver. Hence an advection test on a two-phase flow solver using parallel computations was performed first. The flow domain 
$\varOmega$ is a cube of length 
$1$, and a sphere of radius 
$0.15$ is placed at 
$\boldsymbol {x}=(0.35,0.35,0.35)$. A three-dimensional shear deformation field is defined as
\begin{gather} {U}_1 =2\sin^2({\rm \pi} x)\sin(2{\rm \pi} y)\sin(2{\rm \pi} z)\cos\left(\frac{{\rm \pi} t}{T}\right), \end{gather}
\begin{gather}{U}_2 ={-}\sin^2({\rm \pi} y)\sin(2{\rm \pi} x)\sin(2{\rm \pi} z)\cos\left(\frac{{\rm \pi} t}{T}\right), \end{gather}
\begin{gather}{U}_3 ={-}\sin^2({\rm \pi} z)\sin(2{\rm \pi} x)\sin(2{\rm \pi} y)\cos\left(\frac{{\rm \pi} t}{T}\right), \end{gather}
with time 
$t$ and time period 
$T=3.0$. The 
$\cos ({{\rm \pi} t}/{T})$ term makes the velocity field periodic with respect to time, and ensures that the time integral over a time period at any point in the domain results in zero. This means that any particle moving in the domain will return to its initial position after one time period. Hence it is expected that the sphere will deform under the shear velocity field, get stretched, and then eventually return to its initial shape and position.
Results were obtained using the CLSVOF algorithm on a 
$1/256$ mesh size grid. Figure 1 shows the deformation experienced by the sphere over one time period. The deformed shape at 
$t=1.5$ corresponds to maximum stretching, while at 
$t=3.0$, the sphere has returned to its original position. The CLSVOF algorithm is able to resolve the thin-stretched region at 
$t=1.5$, similar to Ménard, Tanguy & Berlemont (Reference Ménard, Tanguy and Berlemont2007). It can be observed from table 1 that the mass error obtained from our simulations at 
$t=3.0$ when the sphere returns to its original position is comparable with the results reported by Enright et al. (Reference Enright, Fedkiw, Ferziger and Mitchell2002), Ménard et al. (Reference Ménard, Tanguy and Berlemont2007), Wang et al. (Reference Wang, Yang, Koo and Stern2009) and Klitz (Reference Klitz2015).
Figure 1. Evolution of the three-dimensional sphere computed with the CLSVOF method on mesh size 
$1/256$.
Table 1. Mass error of a deformed sphere at 
$t=3.0$.
3.2. Rising bubble in a stratified liquid column
A three-phase flow problem involving a bubble in a stratified liquid column with two liquids having different densities is also used to validate the numerical solver further. An air bubble is placed inside the denser liquid and is allowed to rise gradually and interact with the interface. For an air bubble rising in a stratified liquid column with two liquids, it can either get trapped at the interface of the liquids or penetrate the interface. There is also a possibility that the bubble entrains the heavy-phase liquid if it does penetrate the interface. The condition for a bubble penetrating the interface and the heavier liquid getting entrained, in terms of the radius of the bubble, is given by Greene, Chen & Conlin (Reference Greene, Chen and Conlin1988, Reference Greene, Chen and Conlin1991) and Boyer et al. (Reference Boyer, Lapuerta, Minjeaud, Piar and Quintard2010). Theoretically, the bubble would get trapped at the interface if the radius is 
$r<2.76\ {\rm mm}$, and penetrates the interface if 
$r>2.76\ {\rm mm}$. Simulations of a rising bubble in a stratified column with the physical properties of the fluids the same as those used by Boyer et al. (Reference Boyer, Lapuerta, Minjeaud, Piar and Quintard2010), as mentioned in table 2, are performed for different radii. We have observed that a bubble of radius 
$r=2.00\ {\rm mm}$ gets trapped at the interface, and a bubble of radius 
$r=4.00\ {\rm mm}$ penetrates the interface with very little entrainment of the heavier liquid. For the case with radius 
$r=8\ {\rm mm}$ (see figure 2a), the bubble penetrates the interface while also entraining a large volume of the heavier liquid. This levitated column of heavy liquid elongates and necks down to release a glob of the heavier liquid that is successfully entrained into the lighter liquid, as observed by Greene et al. (Reference Greene, Chen and Conlin1991) in their experimental studies. Figure 2(b) depicts the time evolution of the entrained volume 
$V_e$ of the heavier liquid into a lighter liquid, which is calculated as the amount of the heavy liquid above the initial liquid–liquid interface position. The results obtained match the results of Boyer et al. (Reference Boyer, Lapuerta, Minjeaud, Piar and Quintard2010), which were obtained using the Cahn–Hilliard model.
Table 2. Physical properties of fluids for the rising bubble case.
Figure 2. Simulation of a rising air bubble in a stratified liquid column with radius 
$r=8\ {\rm mm}$. Comparison between present simulations and results obtained by Boyer et al. (Reference Boyer, Lapuerta, Minjeaud, Piar and Quintard2010) for (a) shape of a rising bubble at 
${\rm time} = 0.39$ s, and (b) temporal evolution of entrained volume 
$V_e$ of the heavy liquid in the light liquid.
4. Case set-up
Figure 3(a) shows the computational domain used in the present simulations. A spherical droplet of diameter 
$D$ is placed slightly above the centre of the cubical computational domain of side length 
$4.5D$. The depth of the pool is taken as 
$3D$ in order to ensure that the droplet is sufficiently far away from the computational boundaries. In the experiments (Xie et al. Reference Xie, Forth, Zhu, Helms, Ashby, Shum and Russell2020), the droplets were released from varying heights, but in order to minimize the computational domain size, the droplets are released from a fixed height of 
$0.25D$ but with a different initial velocity. The impact of the droplets is considered for very low Reynolds numbers (
$Re$) and Weber numbers (
$We$). The Reynolds number represents the ratio of inertial forces to viscous forces within a fluid, and the Weber number is the ratio of dynamic pressure (i.e. inertia force) to the surface tension force. The 
$Re$ and 
$We$ values for this investigation are in the ranges 2–43 and 1.8–101, respectively, signifying that splashing and jets are not expected during this impact. All the simulations are performed on grid size 
$1/128$ in all three directions.
Figure 3. (a) Computational domain for the study of hanging droplets. Coacervate model and contact angles at the triple-phase contact line. (b) Interfacial tension at the coacervate with two surfaces and finite coacervate thickness. (c) Representation of the coacervate with a single surface.
The drop liquid is taken as phase 1, the pool liquid is taken as phase 2, and the air is taken as phase 3, for the three-phase flow solver. Figure 3(b) shows the coacervate layer between the drop and pool (immiscible) liquids, where 
$\gamma _{Drop}$ and 
$\gamma _{Pool}$ are the surface tension values for the drop liquid and the pool liquid. The interfacial tensions at the coacervate–drop and coacervate–pool interfaces are given by 
$\gamma _{CD}$ and 
$\gamma _{CP}$. Generally, the coacervate thickness is assumed to be very small, and for the simplicity of modelling, it is taken as a single surface. The two interfacial tensions at the coacervate are combined to give a single interfacial tension at the drop–pool interface (
$\gamma _C=\gamma _{CD}+\gamma _{CP}$), as shown in figure 3(c). The physical properties of the fluids used in the present investigation are similar to those used by Xie et al. (Reference Xie, Forth, Zhu, Helms, Ashby, Shum and Russell2020) and are given in table 3.
Table 3. Physical properties of fluids for the hanging droplet case.
5. Results
5.1. Hanging, intermediate and wrapping droplets
We perform three-dimensional numerical simulations for the above-mentioned configuration and find that the heavier droplet hangs from the lighter liquid interface for certain diameters and impact velocity of the droplet. Figure 4(a) shows the evolution of a hanging droplet (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.137) on the impact of a drop of diameter 
$D=2$ mm, released from height 
$h = 4.98$ mm (impact velocity 
$\sqrt {2gh} = 0.313\ {\rm m}\ {\rm s}^{-1}$). Here, 
$t^*=t/\tau _c$, where capillary time 
$\tau _c$ is defined as 
$\tau _c=\sqrt {\rho _1 D^3/\sigma _{12}}$. As the droplet makes a transition from hanging to sinking, an intermediate case (see supplementary movie 4) is also observed, as shown in figure 4(b) for a droplet of diameter 
$0.6$ mm released from height of 
$h = 61.67$ mm (impact velocity 
$\sqrt {2gh} = 1.1\ {\rm m}\ {\rm s}^{-1}$). Figure 4(c) shows the case for a 
$2$ mm diameter droplet released from height 
$h = 50.97$ mm (impact velocity 
$\sqrt {2gh} = 1\ {\rm m}\ {\rm s}^{-1}$), in which the droplet sinks into the pool upon impact (see supplementary movie 7). It is observed that the droplet begins to slow down even before it makes contact with the pool. As the drop moves closer to the pool, a thin film of air separates the droplet (Duchemin & Josserand Reference Duchemin and Josserand2020) from the pool. This acts as a cushion and is responsible for the decrease in the impact velocity of the droplet. As the droplet makes contact with the pool, the TPCL diameter expands rapidly. After the impact, the droplet loses its kinetic energy drastically by displacing a portion of the pool towards the pool surface. This creates a crater in the pool, shrinking the TPCL diameter. The droplet sits in this crater and hangs from the surface. After the droplet loses all its kinetic energy, it starts moving upwards and keeps oscillating with a very small amplitude until it reaches an equilibrium height. The TPCL diameter again increases during this process. The evolution of the non-dimensional TPCL diameter with respect to the non-dimensional time for different hanging droplet cases is shown in figure 5(a). A capillary wave (Che & Matar Reference Che and Matar2018) is formed upon the impact of the droplet. It is also observed that the equilibrium height and the shape of the hanging droplet are independent of the release height of the droplet as long as it hangs from the surface. This independence of the final shape of the droplet on the impact velocity or release height differs from the results presented by Xie et al. (Reference Xie, Forth, Zhu, Helms, Ashby, Shum and Russell2020) owing to the representation of the coacervate with a single surface.
Figure 4. Evolution of (a) hanging droplet of 2 mm diameter released from 4.98 mm height, (b) intermediate droplet of 0.6 mm diameter released from 61.67 mm height, and (c) wrapping droplet of 2 mm diameter released from 50.97 mm height.
Figure 5. (a) Variation of non-dimensional TPCL diameter with respect to the non-dimensional time for different hanging droplet cases. Here, 
$d^*$ is defined as 
$d^*=d/D$, where 
$d$ is the TPCL diameter. (b) Variation of the non-dimensional height of the droplet for different cases of hanging and wrapping droplets. State diagrams for hanging and wrapping droplets as functions of (c) height and diameter, (d) excess energy and diameter.
Figure 5(b) shows the variation of the height of the centre of mass of the droplet from the pool surface for various cases. It can be observed that there exists a critical depth upon crossing which the droplet gets wrapped. Droplets that do not cross this critical depth tend to hang from the pool surface. The critical depth is found to be 
$0.74$ times the diameter of the droplet. We note that the critical depth is greater than half the diameter of the droplet, i.e. for a brief moment, the droplet goes completely below the pool surface, displacing the pool fluid. Since the computational domain is taken as a closed container such that no fluid exits the domain, the displaced fluid increases the pool height. Increased pool height results in additional pressure head which pushes the droplet upwards. However, if the pool height increases significantly, then it covers the top surface of the droplet and wraps it completely. The droplet sinks when it gets wrapped by the pool fluid.
Simulations for different droplet diameters and release heights are performed (see supplementary movies 1–9). It is observed that the tendency of a droplet to hang from the pool surface increases as the droplet radius or the release height is reduced. A droplet of diameter 
$2$ mm released from height 
$10$ mm gets wrapped and sinks into the pool. In contrast, a droplet of diameter 
$0.4$ mm released from height 
$60$ mm hangs from the pool surface. Figure 5(c) shows both the hanging and wrapping states as functions of the droplet diameter and the release height. A nonlinear curve divides both states. There are also a few cases that lie very close to the curve dividing the two states. These are the cases where a droplet, upon impact with the pool, briefly gets stuck at the pool surface and slowly moves downwards, eventually sinking into the pool. The cases close to the curve dividing the two states are the intermediate cases.
5.2. Energy balance for hanging droplets
Empirical energy calculations are performed to justify hanging and wrapping states for different cases. It is assumed that the droplet is released from height 
$h$, and the entire potential energy is converted into kinetic energy at the time of impact:
\begin{equation} E_{KE}=\frac{\rm \pi}{6}\,D^3\rho_1 gh. \end{equation}
This is the entire energy available with the droplet that is used to overcome different forms of energy requirements. Three different forms of energy loss are considered here for energy balance. First, a part of the available energy is spent to displace the pool fluid to the pool surface to create a crater for the droplet. It is observed from figure 5(b) that if a droplet is getting wrapped, then it needs to attain a critical depth. From this, the displaced volume 
$\Delta \mathcal {V}$ is approximated as
\begin{equation} \Delta \mathcal{V} \approx \frac{\rm \pi}{4}\,D^2\times 0.74D + \frac{\rm \pi}{12}\,D^3. \end{equation}Here, it is assumed that the crater is cylindrically shaped with a hemisphere at one of its ends. The centre of mass of the crater is given as
\begin{equation} z_c = \frac{\dfrac{\rm \pi}{4}\,D^2\times 0.74D \times \dfrac{0.74}{2}\,D + \dfrac{\rm \pi}{12}\,D^3\times \left(0.74+\dfrac{3}{16}\right)D}{\dfrac{\rm \pi}{4}\,D^2\times 0.74D + \dfrac{\rm \pi}{12}\,D^3}. \end{equation}
This 
$z_c = 0.5431D$ is the height by which the crater needs to be lifted, and the energy required for this is calculated as potential energy loss. Taking a correction factor 
$k_{PE}$, the potential energy loss is evaluated as
\begin{equation} E_{PE} \approx k_{PE} \times \rho_2\,\Delta \mathcal{V}\,gz_c = k_{PE} \times 4.4916\rho_2D^4. \end{equation}The correction factor for the potential energy loss is taken as unity. By taking the shape of the crater as defined above, change in surface area can be evaluated for different surfaces. Taking the product of these surface changes with their respective surface tension values gives an estimate of the energy required for the destruction and creation of new surfaces:
\begin{equation} E_{S} \approx \left( 0.74{\rm \pi} D^2\sigma_{23} + \frac{\rm \pi}{2}\,D^2\sigma_{12}\right)-\left(\frac{\rm \pi}{2}\,D^2\sigma_{13} + \frac{\rm \pi}{4}\,D^2\sigma_{23} \right). \end{equation}
Using the values of 
$\sigma _{12}$, 
$\sigma _{23}$ and 
$\sigma _{13}$, this energy is further approximated as
\begin{equation} E_{S} \approx k_S \times 0.005{\rm \pi} D^2. \end{equation}
The actual crater is not exactly cylindrical, but rather has curved edges and capillary waves. Therefore, the actual surface generated should be bigger than estimated. As a result, the adjustment factor 
$k_S$ should be greater than 
$1$. The liquid in the pool is highly viscous, hence large viscous losses are also expected due to the motion of the droplet into the pool. Since the force experienced by a droplet when moving through another fluid is not known exactly, the following assumptions are made to approximate this energy loss: (a) the droplet is assumed to be a rigid sphere moving through the pool; (b) flow speed past the droplet is taken as constant. Under these assumptions, the drag force experienced by the droplet is given by
\begin{equation} F_d = \tfrac{1}{2}C_d\rho_2 V^2 \times A, \end{equation}
where 
$A={\rm \pi} D^2/4$ is the frontal area. The distance travelled by the droplet 
$z_d$ is the sum of the height of the centre of mass of the drop at the time of impact above the interface (
$0.5D$) and the critical depth (
$0.74D$). We also have to include a correction factor 
$k_V$ to accommodate the aforementioned assumptions. This gives the viscous losses as
\begin{equation} E_V \approx k_V \times \tfrac{1}{2}C_d\rho_2 V^2 A z_d = 0.155 k_V C_d\rho_2 V^2{\rm \pi} D^3. \end{equation}
When the droplet descends from the interface, its velocity decreases, and its shape changes, resulting in a decrease in the drag coefficient. Therefore, the overall viscous losses will be lower than the estimated value and should also be accounted for by the correction factor. Here, 
$C_d$ is the drag coefficient for flow past a sphere and is dependent on the droplet diameter 
$D$, impact velocity 
$V$, density 
$\rho _2$, and viscosity 
$\mu _2$ of the pool. The Reynolds number 
$Re = {\rho _2 V D}/{\mu _2}$ satisfies 
$1< Re<1000$. Hence the drag coefficient is defined using the relation given by Schiller and Naumann (Flemmer & Banks Reference Flemmer and Banks1986). Considering the above-mentioned four energies, it is determined whether a droplet, if it crosses the critical depth, still has additional energy to move further downwards. Excess energy is calculated as
\begin{equation} E_{{excess}} = E_{KE}-E_{PE}-E_S-E_V. \end{equation}
The values of the correction factors 
$k_S$ and 
$k_V$ are tuned such that the available data set for the final state of droplet impact gives a distinct distribution in terms of the excess energy. Figure 5(d) shows the state diagram for hanging and wrapping droplets as functions of excess energy and the diameter of the droplet. It can be seen that the droplets with sufficient energy to spend on different losses tend to get wrapped and sink into the pool, whereas the droplets that have less energy to begin with, such that they have negative excess energy, tend to hang from the pool surface. There are also intermediate cases where the available energy is nearly equal to the energy required and thus has close to zero excess energy. These droplets initially lose their entire kinetic energy upon impact and then gradually sink into the pool. It is observed that the majority portion of the available energy is spent to overcome the viscous loss, and the remaining energy is spent for surface energy. A very small part of the available energy is spent on potential energy loss. Thus a larger droplet with higher initial energy can still hang from the surface if either the viscosity of the pool fluid is increased or the interfacial tension value used for the coacervate is increased. We note here that the excess energy is just a function of 
$D$ and 
$h$, and it converts the nonlinear distribution of hanging and wrapping droplets in figure 5(c) into a linear distribution in figure 5(d).
5.3. Force balance for hanging droplets
The viscous force has a significant impact on the droplet's rate of descent; nevertheless, after it has reached equilibrium, it is the surface tension force and the buoyant forces that are responsible for maintaining the droplet's attachment to the surface by balancing its weight. We present the calculations for the force balance based on an analytical approach and the outcomes of the numerical simulations. The weight of the droplet is calculated as
\begin{equation} F_{w} = \mathcal{V}_1 \times \rho_1 \times g. \end{equation}
Here, 
$\mathcal {V}_1$ is the total volume of the droplet. The buoyant force is defined as
\begin{equation} F_{b} = \mathcal{V}_2 \times \rho_2 \times g. \end{equation}
Here, 
$\mathcal {V}_2$ is the volume of pool fluid displaced by the droplet below the TPCL. It is worth noting that the droplet at equilibrium is not completely immersed beneath the pool's surface. A little portion of the droplet remains above the TPCL line. There are also pockets of air bubbles trapped between the droplet and pool interfaces. The volume of these air bubbles is also included in 
$\mathcal {V}_2$ to calculate the buoyant forces. Now consider a system with a droplet including the droplet–air interface and the droplet–pool interface. The forces acting on the system in the vertical direction are only the gravitational force, the buoyant force, and the surface tension due to the air–drop interface and drop–pool interface. Since the hanging droplet system consists of two different surfaces wrapped around a common ring, i.e. the TPCL, the surface tension forces due to each of the two interfaces can be evaluated by taking the product of the pressure jump across the interface and the projected area at their boundary. The pressure jump at the interface can be calculated using Young's Laplace equation:
\begin{equation} \Delta p_{ij} = \frac{2\sigma_{ij}}{R_{ij}}, \end{equation}
where 
$R_{ij}$ is the radius of curvature of the interface between the phases 
$i$ and 
$j$. Therefore, the vertical force on the hanging droplet owing to the surface tension computed from the values of 
$R_{12}$ and 
$R_{13}$ obtained from the simulations is given by
\begin{equation} F_{\sigma} = \left(\frac{2\sigma_{12}}{R_{12}}-\frac{2\sigma_{13}}{R_{13}}\right) \times \frac{{\rm \pi} d^2}{4}. \end{equation}
Based on the values of 
$\mathcal {V}_1$ and 
$\mathcal {V}_2$ obtained from our simulations, in order to satisfy the force balance on the hanging droplet in the vertical direction, the ideal value of the vertical component of surface tension force should be
\begin{equation} F_{\sigma0} = F_w-F_b. \end{equation}
The 
$F_{\sigma }$ values computed from the simulations are in fact very close to 
$F_{\sigma 0}$, as demonstrated in table 4 for four cases of hanging droplets signifying the dynamical balance.
Table 4. Calculation for surface tension force balance for the hanging droplet.
Our simulations show that the shape of a hanging droplet at equilibrium resembles a combination of two spherical caps with varying radii. Hence we model the droplet with two spherical sections of radii 
$R_{12}$ and 
$R_{13}$, respectively, as shown in figure 6(a). Considering the two interfaces as part of purely spherical sections, the radius of curvature and the TPCL diameter can be related as
\begin{equation} d = 2R_{12}\sin\beta = 2R_{13}\sin\phi. \end{equation}
Using (5.15), the volumes of the upper and lower spherical sections, 
$V_{13}$ and 
$V_{12}$, can be evaluated as
\begin{gather} V_{13} = \frac{\rm \pi}{24\sin^3\phi}\,d^3(2+\cos\phi)(1-\cos\phi)^2,\end{gather}
\begin{gather}V_{12} = \frac{\rm \pi}{24\sin^3\beta}\,d^3(2-\cos\beta)(1+\cos\beta)^2. \end{gather}The droplet and the pool fluids are taken as immiscible because there is no chemical reaction taking place at the interface. Therefore, the volume of the droplet must be conserved:
\begin{equation} V_{12}+V_{13} = \frac{\rm \pi}{6}\,D^3. \end{equation}Again using (5.15) in (5.13), the vertical surface tension force on the droplet is calculated analytically as
\begin{equation} F_{\sigma,model} = {\rm \pi}d (\sigma_{12}\sin\beta-\sigma_{13}\sin\phi). \end{equation}However, considering the vertical force balance on the droplet, i.e. using (5.14), (5.17) and (5.16), we get
\begin{equation} F_{\sigma,model} = \left((\rho_1-\rho_2)V_{12}+(\rho_1-\rho_3)V_{13}\right)g. \end{equation}Apart from the force balance on the droplet, the interfaces between the three phases, air, droplet and pool, are also considered to be massless. Therefore, at the junction of the three phases, i.e. at the TPCL, the vertical and the horizontal surface tension forces due to the three interfaces must balance each other (Phan et al. Reference Phan, Allen, Peters, Le and Tade2012; Xie et al. Reference Xie, Forth, Zhu, Helms, Ashby, Shum and Russell2020). Therefore,
\begin{gather} \sigma_{23}\sin\alpha-\sigma_{12}\sin\beta+\sigma_{13}\sin\phi = 0, \end{gather}
\begin{gather}\sigma_{23}\cos\alpha+\sigma_{12}\cos\beta-\sigma_{13}\cos\phi = 0. \end{gather}
Eliminating 
$\alpha$ (see figure 3b) from (5.21) and (5.22), we get
\begin{equation} (\sigma_{13}\sin\phi-\sigma_{12}\sin\beta)^2 + (\sigma_{13}\cos\phi-\sigma_{12}\cos\beta)^2 = \sigma^2_{23}. \end{equation}
Equations (5.18), (5.19), (5.20) and (5.23) can be solved simultaneously to obtain the values for 
$\beta$, 
$\phi$, 
$d$ and 
$F_{\sigma,model}$. Figures 6(b), 6(c) and 6(d) show an excellent match of the contact angles (
$\beta$ and 
$\phi$), radii of curvature for the two interfaces, and the surface tension forces respectively between our simulations and the analytical solution obtained by solving the force balance equations for the hanging drops.
Figure 6. (a) Shape of droplet at equilibrium. (b) Contact angles 
$\beta$ and 
$\phi$ versus diameter of the droplet at equilibrium. (c) Radius of curvature for the interface versus diameter of the droplet at equilibrium. (d) Surface tension force on the droplet versus diameter of the droplet at equilibrium.
6. Conclusions
In this work, the impact of a droplet into a pool of immiscible liquid is investigated using three-dimensional three-phase flow simulations. The results from the numerical simulations suggest that the droplet upon impact can either hang from the liquid surface or get wrapped into the pool and sink eventually. In some rare cases, the droplet even gets stuck at the interface and gradually sinks into the pool. All three cases obtained from the numerical results are shown to happen in experiments of Xie et al. (Reference Xie, Forth, Zhu, Helms, Ashby, Shum and Russell2020) as well. Further, a parametric study of the droplet impact is done to understand the effect of droplet diameter and release height on the final state of the droplet. It is observed that a nonlinear curve in terms of droplet diameter and release height separates the hanging and wrapping states. As the droplet diameter or the release height is increased, the droplets move from the hanging state to the wrapping state. It is observed that the shape of the droplet at equilibrium does not vary with release height for hanging droplets. A hanging droplet of a given diameter tends to have a unique final state. This behaviour of the hanging droplets is different from the observation of Xie et al. (Reference Xie, Forth, Zhu, Helms, Ashby, Shum and Russell2020). The simplicity of the model used for the coacervate is the probable reason for this divergence from the experimental results. This suggests that the coacervate needs more sophisticated modelling for its physical properties even if it is considered to have no mass.
Further, an approximate energy balance is presented for the droplet impact. It is shown that the major portion of the kinetic energy available with the droplet is dissipated as viscous losses. The remaining energy is converted to the surface and potential energy owing to the formation of the crater. The energy balance is then used to determine whether a given heavier droplet will float or sink in the pool of the lighter liquid. Additionally, we solve the force balance equations of a hanging drop analytically at the equilibrium position. The values of the contact angle, the radii of curvature, and the force due to surface tension at the TPCL obtained from this dynamical balance show an excellent match with the simulations.
A natural extension of this work will be to perform a parametric study to understand the effect of other fluid parameters such as the viscosity and the surface tension values. Furthermore, larger-sized droplets with much higher energy can be simulated to potentially get some new states like droplets breaking off the surface leaving a secondary droplet at the surface.
Supplementary movies
Supplementary movies of hanging, intermediate and wrapping droplets are available at https://doi.org/10.1017/jfm.2023.137.
Acknowledgements
The support and the resources provided by PARAM Sanganak under the National Supercomputing Mission of the Government of India at the Indian Institute of Technology, Kanpur are gratefully acknowledged.
Funding
We gratefully acknowledge the support of the Science and Engineering Research Board, Government of India grant no. SERB/ME/2020318. We also want to thank the Office of Research and Development of the Indian Institute of Technology, Kanpur for the financial support through grant no. IITK/ME/2019194.
Declaration of interests
The authors report no conflict of interest.
Author contributions
The authors contributed equally to analysing data and reaching conclusions, and in writing the paper.
