1. Introduction
Viscoelastic fluids are a class of non-Newtonian fluids that exhibit a viscous and an elastic response to the application of stress (Snoeijer et al. Reference Snoeijer, Pandey, Herrada and Eggers2020). Some examples include biological fluids such as saliva (Mitchinson Reference Mitchinson2010), the children’s toy Silly Putty and other complex fluids such as surfactant solutions and emulsions. The common feature that is present in these different viscoelastic fluids is that they all contain molecules that are able to stretch and relax over some time scale. Polymer solutions are a common viscoelastic fluid and have been used extensively to study the influence of material properties on flow behaviour and the many interesting phenomena that arise as they flow (Datta et al. Reference Datta2022). For example, the addition of polymers to a solution leads to drag reduction in turbulent pipe flows (Toms Reference Toms1948; Berman Reference Berman1978) and can alter drop breakup by producing long and stable threads (Anna & McKinley Reference Anna and McKinley2001).
More recently, the coalescence of drops of polymer solution has been studied under different configurations varying both the drop geometry and the substrate wettability (Varma et al. Reference Varma, Saha, Mukherjee, Bandopadhyay, Kumar and Chakraborty2020, Reference Varma, Saha and Kumar2021; Chen et al. Reference Chen, Pan, Nie, Ma, Fang and Yin2022; Dekker et al. Reference Dekker, Hack, Tewes, Datt, Bouillant and Snoeijer2022; Varma et al. Reference Varma, Dasgupta and Kumar2022a
, Reference Varma, Rajput and Kumar2022b
; Sivasankar et al. Reference Sivasankar, Etha, Hines and Das2023; Eggers, Sprittles & Snoeijer Reference Eggers, Sprittles and Snoeijer2024; Rostami et al. Reference Rostami, Erb, Azizmalayeri, Steinmann, Stark and Auernhammer2025). Here, we focus on the coalescence of sessile drops (figure 1). Studies have shown that in the case of semi-dilute polymer solutions and inertially dominated coalescence, which often correspond to a large contact angle
$\theta$
, the temporal evolution of the bridge height
$h_0$
(see figure 1) follows the same scaling as its Newtonian counterpart:
$h_0 \propto t^{2/3}$
(Varma et al. Reference Varma, Dasgupta and Kumar2022a
; Dekker et al. Reference Dekker, Hack, Tewes, Datt, Bouillant and Snoeijer2022). While the temporal scaling of the bridge height was unaffected, the shape of the interface was highly altered by the polymer stress. The Newtonian, self-similar scaling that describes the evolution of the interface at early times after coalescence breaks down for these semi-dilute polymer solutions and different self-similar scalings were proposed taking into account the stress induced by the polymers (Dekker et al. Reference Dekker, Hack, Tewes, Datt, Bouillant and Snoeijer2022; Varma et al. Reference Varma, Dasgupta and Kumar2022a
). As the polymer concentration increased, the magnitude of the temporal scaling exponent decreased from
${2}/{3}$
to about
$0.5$
in certain cases and even lower in others (Varma, Saha & Kumar Reference Varma, Saha and Kumar2021, Reference Varma, Dasgupta and Kumar2022a
; Chen et al. Reference Chen, Pan, Nie, Ma, Fang and Yin2022; Rostami et al. Reference Rostami, Erb, Azizmalayeri, Steinmann, Stark and Auernhammer2025). The nature of the change in the scaling exponent was shown to be affected by the way in which the drops were deposited. When two drops of fixed volume are deposited by allowing them to impact a substrate and spread until they come into contact, the scaling exponent decreases to approximately
$0.5$
with increasing polymer concentration. When fluid is injected into two drops so that their volume increases until they touch and coalesce, the scaling exponent becomes even lower (Varma et al. Reference Varma, Dasgupta and Kumar2022a
). While much is known about the effect of elasticity on the coalescence dynamics, a comprehensive understanding of the departure of these dynamics from Newtonian dynamics is still lacking. To the best of our knowledge, the coalescence of viscoelastic sessile drops in the small contact angle regime also remains unexplored.
The shape of the interface during the coalescence of Newtonian and polymeric drops. (a) Schematic of the side-view profile of two drops during a typical coalescence experiment. Experimental images of large-
$\theta$
coalescence of (b) water and (c) 0.5 wt% PEO drops show a significant difference in the shape of the interface. On the other hand, images of small-
$\theta$
coalescence of (d) 1000 cSt silicone oil and (e) 0.5 wt% PEO drops show similar shape of the interface. Scale bars represent 0.1 mm.
Here, we combine laboratory experiments and numerical simulations using the Oldroyd-B model to investigate the coalescence of sessile viscoelastic drops at small and large static advancing contact angle of the drops at the moment of coalescence, hereafter referred to simply as the contact angle
$\theta$
. Experimentally, we consider the coalescence of sessile semi-dilute polymeric drops (polyethylene oxide (PEO); molecular weight
$M_w=4 \times 10^6$
g mol–1). To illustrate the influence of the contact angle, figure 1 shows side-view profiles of the coalescence of drops with and without polymer at large and small contact angles. The difference in the shape of the interface with and without polymer is clearly evident at large contact angles, as shown in figure 1(b,c), and as previously reported (Dekker et al. Reference Dekker, Hack, Tewes, Datt, Bouillant and Snoeijer2022). In contrast, the shape of the interface with and without polymer at small contact angles looks nearly the same, as shown in figure 1(d,e).
In what follows, in § 3, we analytically and experimentally show that the coalescence dynamics is virtually unaffected by the presence of polymers in the limit
$\theta \ll 1$
, which we show is connected to the Deborah number,
$\textit{De}$
, being small. The Deborah number is the ratio between the relaxation time of the polymer and the characteristic time scale of the background flow. In § 4, we use numerical simulations to study coalescence at
$\theta \gtrapprox 1$
. We show that the Oldroyd-B model can capture the effect of polymers seen in the experiments and further illustrate the coalescence dynamics as a function of the Deborah number and the elastocapillary number,
$\textit{Ec}$
, of the polymer solution. The elastocapillary number is the ratio between the shear modulus of the polymer solution and the characteristic stress in the fluid. Exploring the various limits of the Oldroyd-B model, we explain our numerical results and provide more insight into the departure of the coalescence dynamics from the Newtonian case. In both our experiments and simulations, we make sure that the time scale of coalescence is clearly separated from all other relevant time scales. For instance, in the experiments, the droplets spread slowly over several minutes prior to contact, while the actual coalescence occurs over just a few seconds. In simulations, this feature is straightforward to achieve by initialising all velocities to zero.
2. Governing equations
Here, we introduce the general equations that will be rescaled later to identify the dimensionless parameters for the coalescence problem. The governing mass and momentum equations for incompressible flows are given by
\begin{gather} \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot }\tilde {\boldsymbol{u}}=0, \end{gather}
\begin{gather} \rho \frac {\textrm {D} \tilde {\boldsymbol{u}}}{\textrm {D} \tilde {t}} = -\tilde {\boldsymbol{\nabla }} \tilde {p} + \tilde {\boldsymbol{\nabla }} \boldsymbol{\cdot }\tilde {\boldsymbol{\tau }}, \end{gather}
where
${\textrm {D} }/{\textrm {D} \tilde {t}}$
denotes the material derivative; the tilde is used to represent dimensional variables. Here,
$\tilde {\boldsymbol{u}}$
is the velocity field,
$\tilde {p}$
is the pressure,
$\rho$
is the density and
$\tilde {\boldsymbol{\tau }}$
is the deviatoric stress tensor. For polymeric fluids, the deviatoric stress tensor can be written as the sum of the solvent and polymeric contributions:
$\tilde {\boldsymbol{\tau }} = \tilde {\boldsymbol{\tau }}^s + \tilde {\boldsymbol{\tau }}^{p}$
. The solvent contribution is the Newtonian stress
$\tilde {\boldsymbol{\tau }}^s = 2\mu_s \tilde {\boldsymbol{E}}$
, where
$\mu_s$
is the solvent viscosity and
$\tilde {\boldsymbol{E}} = ({1}/{2})[(\tilde {\boldsymbol{\nabla }} \tilde {\boldsymbol{u}})+(\tilde {\boldsymbol{\nabla }} \tilde {\boldsymbol{u}})^T]$
is the rate of strain tensor. The polymer contribution is modelled by the Oldroyd-B model according to
\begin{gather} \tilde {\boldsymbol{\tau }}^{p} = G (\boldsymbol{A} - \boldsymbol{I}), \end{gather}
\begin{gather} \overset {\mathtt {\boldsymbol{\nabla }}}{\boldsymbol{A}} = -\frac {1}{\lambda }(\boldsymbol{A} - \boldsymbol{I}), \end{gather}
\begin{gather} \lambda \overset {\mathtt {\boldsymbol{\nabla }}}{\tilde {\boldsymbol{\tau }}^p} + \tilde {\boldsymbol{\tau }}^p = 2 G \lambda \tilde {\boldsymbol{E}}, \end{gather}
where
$G$
is the elastic modulus of the polymer solution,
$\boldsymbol{A}$
is the conformation tensor of the polymer,
$\boldsymbol{I}$
is the identity tensor,
$\lambda$
is the relaxation time of the polymer solution and
$G \lambda$
is the polymeric viscosity. The upper convected derivative is defined as
$\overset {\mathtt {\boldsymbol{\nabla }}}{\boldsymbol{A}} = (\mathrm{D}/{{\mathrm{D}} \tilde {t}})\boldsymbol{A} - (\tilde {\boldsymbol{\nabla }} \tilde {\boldsymbol{u}})^T \boldsymbol{\cdot }\boldsymbol{A} - \boldsymbol{A} \boldsymbol{\cdot }\tilde {\boldsymbol{\nabla }} \tilde {\boldsymbol{u}}$
. Note that (2.2a
) is a state of stress assumption and (2.2b
) is an evolution equation for the microstructure of the polymer chains, which are allowed to stretch/compress under flow and relax with the longest relaxation time
$\lambda$
(Stone, Shelley & Boyko Reference Stone, Shelley and Boyko2023). Equation (2.2c
) is usually referred to as the Oldroyd-B model, and can be obtained by taking the upper convected derivative of (2.2a
) and appropriately substituting (2.2b
) (Snoeijer et al. Reference Snoeijer, Pandey, Herrada and Eggers2020).
3. The small contact angle limit
We start by considering the coalescence of polymeric drops in the small contact angle or thin-film limit. We identify the important dimensionless parameters governing the dynamics and show experimental results from imaging the three-dimensional shape of the interface during coalescence.
3.1. Dimensionless equations and parameters
In the small contact angle limit,
$\theta \ll 1$
, the lubrication approximation has been shown to successfully describe the dynamics of coalescing sessile Newtonian drops (Ristenpart et al. Reference Ristenpart, McCalla, Roy and Stone2006; Hernández-Sánchez et al. Reference Hernández-Sánchez, Lubbers, Eddi and Snoeijer2012; Kaneelil et al. Reference Kaneelil, Pahlavan, Xue and Stone2022). Therefore, the equations can be rescaled using the following variables:
$x=\tilde {x}/\ell _c,\,y=\tilde {y}/\ell _c,\,z=\tilde {z}/(\theta \ell _c),\,u_x=\tilde {u_x}/u_c,\,u_y=\tilde {u_y}/u_c,\,u_z=\tilde {u_z}/(\theta u_c),\,t=\tilde {t}/(\ell _c/u_c)$
and
$p=\tilde {p}\theta ^2 \ell _c/(\mu u_c)$
, where
$x$
and
$y$
are the in-plane coordinates on the substrate and
$z$
is the out-of-plane coordinate, as shown in figure 2(a). Also,
$\ell _c=\sqrt {\gamma /(\rho g)}$
is the capillary length, which is a good approximation for the largest height of the drop in our experiments,
$\mu$
is the viscosity of the solution, and
$u_c=\gamma \theta ^3/(3\mu )$
is the velocity scale for coalescence (Hernández-Sánchez et al. Reference Hernández-Sánchez, Lubbers, Eddi and Snoeijer2012). Also,
$\textit{Oh}_{\theta }=\mu /\sqrt {\rho \gamma \ell _c}$
denotes the Ohnesorge number. Note that we neglect inertial effects in this limit, as usual, since the inertial terms in the momentum equation (2.1b
) will vary, according to the above scalings, as
$\theta ^5/\textit{Oh}_{\theta }^2$
(about
$10^{-6}$
in our experiments) and will always be small when
$\theta \ll 1$
.
Three-dimensional reconstruction of the shape of the interface using FS-SS imaging. (a) Schematic showing the experimental set-up and the drop geometry. (b) Sequence of experimental images showing a reference frame taken before the drop appeared, and two time steps during the spreading and coalescing of 1 wt% PEO drops. (c) The three-dimensional reconstruction of the interface shape corresponding to
$t=1$
s after coalescence.
Rescaling the Oldroyd-B equations with these variables gives
\begin{gather} \boldsymbol{\tau }^p=3 \textit{Ec}_{\theta }(\boldsymbol{A}-\boldsymbol{I}),\end{gather}
\begin{gather} \overset {\mathtt {\boldsymbol{\nabla }}}{\boldsymbol{A}}=-\frac {1}{De_\theta }(\boldsymbol{A}-\boldsymbol{I}), \end{gather}
where the Deborah number and the elastocapillary number are defined as
\begin{equation} De_{\theta } = \frac {\lambda \gamma \theta ^3}{3\mu \ell _c} \quad\textrm {and}\quad \textit{Ec}_{\theta }=\frac {G\ell _c}{\theta \gamma }, \end{equation}
with the subscript
$\theta$
indicating the
$\theta \ll 1$
limit. Similarly, we rescale (2.2c
) and focus on the
$z$
components of the equations, which become
\begin{gather} De_\theta \,\overset {\mathtt {\boldsymbol{\nabla }}}{\tau ^p_{zj}} + \tau ^p_{zj} = 3\, \frac {\textit{Ec}_{\theta } De_\theta }{\theta }\, \frac {\partial u_j}{\partial z} + {\mathcal O}(\theta ^2),\quad j = x,y, \end{gather}
\begin{gather} De_\theta \,\overset {\mathtt {\boldsymbol{\nabla }}}{\tau ^p_{zj}} + \tau ^p_{zj} = 6\, \textit{Ec}_{\theta } De_\theta \, \frac {\partial u_j}{\partial z},\quad j = z. \end{gather}
Thus, the dimensionless parameters that affect the dynamics of the interface can be identified as
$\theta ,\ De_\theta $
and
$\,\textit{Ec}_{\theta }$
. In these shear-dominated thin-film flows, the largest velocity gradients and therefore the largest stresses are in the
$xz$
and
$yz$
directions, corresponding to (3.3a
). From the definition of the dimensionless parameters and (3.3a
), we can observe for
$\theta \ll 1$
that
\begin{equation} \tau ^p_{zj} \sim \theta \,\frac {\partial u_j}{\partial z} +{\mathcal O}(\theta ^2), \quad j=x,y , \end{equation}
which suggests that polymer stress will not affect coalescence dynamics at leading order. The experiments described below agree well with this prediction. The result can be attributed to the significant dependence of the Deborah number on the contact angle, i.e.
$De_\theta \sim \theta ^3$
. This relationship results from the characteristic velocity scaling similarly,
$u_c \sim \theta ^3$
, which is a feature of capillary-driven thin-film flows. The same result can be immediately seen by taking the
$\theta \rightarrow 0$
limit of (3.2), where
$Ec$
becomes infinite and
$De$
goes to zero. This is the Newtonian limit where there is no memory and the stress response is instantaneous.
3.2. Experimental set-up
3.2.1. Materials and methods
The experiments were performed using PEO (Sigma-Aldrich; molecular weight
$M_w=4 \times 10^6\,\mathrm{g\,mol}^{-1}$
) dissolved in deionised water at three different concentrations:
$c= 0.1$
,
$0.5$
and
$1.0$
wt%. In terms of the critical coil overlap concentration
$c^*=1/(0.072 \,M_w^{0.65})$
above which the polymer coils may overlap one another (Tirtaatmadja, McKinley & Cooper-White Reference Tirtaatmadja, McKinley and Cooper-White2006), these solutions correspond to
$c/c^*=1.4,\ 7.0$
and
$14.1$
, respectively. Note that the units of
$c^*$
are g ml−1 and therefore the prefactor of the power law has units (ml mol
$^{0.65}$
) (g
$^{-1.65}$
). After the solutions were made, they were mixed with a magnetic stirrer for at least 24 h. The surface tensions of the three solutions were measured using the pendant drop method and were
$\gamma = [62.7 \pm 2.5,\,62.8 \pm 1.9,\,63.1 \pm 2.7]$
mN m−1, respectively. The shear viscosities
$\mu$
of the three solutions, at low shear rates, were measured using a rheometer (Anton Parr MCR 302e) and were
$\mu = [2.37 \pm 0.24,\,27.5 \pm 0.7, 263.4 \pm 8.8]$
mPa s. The relaxation times for the solutions were obtained from the literature (Dekker et al. Reference Dekker, Hack, Tewes, Datt, Bouillant and Snoeijer2022), where they were measured from the thinning dynamics of the neck during pinch-off of a PEO drop.
We used glass microscope slides as the substrate for the experiments. They were cleaned by sequentially immersing and sonicating the slides for 15 min each in a surfactant solution, deionised water, ethanol and acetone bath. A clean substrate is very important to prevent pinning of the triple line and to make sure that the two drops have the same contact angle upon contact and coalescence. The contact angle
$\theta$
of the drops with the substrate at the moment of coalescence was within
$7.5^{\circ }$
–
$12^{\circ }$
(
$0.13$
–
$0.21\,\textrm {rad}$
). The drops were dispensed successively through a needle and were placed far enough apart to spread and reach small angles before contact. The angles of the drops were measured a posteriori from the experimental images to confirm symmetric coalescence.
3.2.2. Imaging and processing
In order to image the three-dimensional shape of the interface during coalescence, we used free-surface synthetic schlieren (FS-SS) imaging (Moisy, Rabaud & Salsac Reference Moisy, Rabaud and Salsac2009). In this set-up, a dot pattern is placed underneath the substrate that holds the drops and the imaging is done with a camera from the top (Phantom V7.3; 400 fps). Figure 2(a) shows a schematic of the FS-SS set-up. This technique utilises the fact that the curved shape of the drop acts as a lens to distort the image of the dot pattern that is under the drop. By measuring the distortion, or the displacement field of the dot pattern compared with a reference image, we can calculate the shape of the interface that caused the distortion. Figure 2(b) shows a reference image and two images at
$t=0$
and
$t=3$
s during the coalescence of 1 wt% PEO drops. Notice that the latter two images show a distorted dot pattern relative to the reference image. Figure 2(c) shows the reconstructed three-dimensional shape of the interface from the experimental images. More information about how we calculate the interface profile from the images can be found in Appendix D.
3.3. Results and discussion
3.3.1. Height of the bridge
In coalescence problems, analysing the height of the bridge is a starting point for understanding the dynamics of the system (Hernández-Sánchez et al. Reference Hernández-Sánchez, Lubbers, Eddi and Snoeijer2012; Kaneelil et al. Reference Kaneelil, Pahlavan, Xue and Stone2022). Note that the height of the bridge,
$h_0(t)$
, refers to the height of the interface at the point where the two drops initially made contact. Figure 3(a) summarises the time evolution of the height of the bridge for different polymer concentrations. Data from multiple experiments are shown for each polymer concentration and span the ranges
$De_\theta = [0.002,\,0.06]$
and
$\textit{Ec}_\theta = [0.07,\,1.7]$
. Theory for the evolution of the height of the bridge for viscous Newtonian drops predicts a linear scaling with time,
$h_0(t) =vt= (A\theta u_c) t=A ({\gamma \theta ^4}/{3 \mu }) t$
, where
$v$
is the velocity in the
$z$
direction and
$A \approx 0.818$
is a prefactor that can be determined uniquely (Hernández-Sánchez et al. Reference Hernández-Sánchez, Lubbers, Eddi and Snoeijer2012). Note that the variability in the
$h_0(t)$
data for a given polymer concentration in figure 3(a) arises from the difference in
$\theta$
between experiments, which leads to varying vertical velocities of the interface. We fitted the data using a power law,
$h_0(t) \propto t^\alpha$
, and the results are shown in figure 3(b). The average power-law exponents
$\alpha$
for the 0.1, 0.5 and 1 wt% PEO cases were
$1.14\pm 0.04,\, 1.16\pm 0.05$
and
$1.03\pm 0.08$
, respectively. The coefficients of determination in all cases were above
$R^2=0.99$
.
The time evolution of the height of the interface
$h_0(t)$
at the initial coalescence point. (a) Raw data showing
$h_0(t)$
from experiments using three different polymer concentrations, spanning the ranges
$De_\theta = [0.002,\,0.06]$
and
$\textit{Ec}_\theta = [0.07,\,1.7]$
. (b) Average power-law exponent
$\alpha$
from fitting the data for the different polymer concentrations. (c) The
$h_0$
versus
$t$
data rescaled according to Newtonian viscous scaling. Rescaling reasonably collapses the data, and the black line has a power-law exponent
$\alpha =1$
and a prefactor
$A=0.818$
, predicted by the viscous theory.
We rescaled the height evolution data according to Newtonian scaling, as shown in figure 3(c), and observed a reasonable collapse of the data. The rescaled data are expected to have a slope of
$A$
, which is the slope of the black line plotted in figure 3(c). Although the data show reasonable collapse compared with raw data, there exists some variability. We believe that this variability originates from the error in experimental measurement of the contact angle of the drops, which gets amplified due to the
$\theta ^4$
dependence. We now reveal the shape of the interface and show that it is also in good agreement with the Newtonian counterpart of this experiment.
3.3.2. Shape of the interface
During coalescence of sessile drops, the shape of the interface that forms near the coalescence point is three-dimensional and resembles a saddle (Kaneelil et al. Reference Kaneelil, Pahlavan, Xue and Stone2022). Unlike the case of coalescence of spherical drops, the three-dimensionality arises from the presence of the contact line. In the
$xz$
plane, the height at the coalescence point
$h_0(t)$
will be the lowest point on the interface and the interface will slope upwards to form the edge of the drop. In the
$yz$
plane,
$h_0(t)$
will be the highest point on the interface and the interface will slope downwards on either side towards the contact points on the substrate (see figures 2
a and 4
a,d). Figure 4(b) shows the dynamic shape of the interface in the
$xz$
plane,
$h(x,y=0,t)$
, as two drops of 0.5 wt% PEO coalesce. In the Newtonian case, the height profile in the
$xz$
plane has a self-similar profile where both
$h(x,y=0,t)$
and
$x$
are rescaled by
$h_0(t)$
(Hernández-Sánchez et al. Reference Hernández-Sánchez, Lubbers, Eddi and Snoeijer2012). The results in figure 4(c) show that this self-similar scaling collapses the data well.
The interface profiles along the
$x$
and
$y$
axes from the coalescence of 0.5 wt% PEO drops with
$\theta \approx 8.1^{\circ }$
, corresponding to
$De_\theta = 0.009$
and
$\textit{Ec}_\theta = 0.74$
. (a) Schematic of the interface in the
$xz$
plane where the height
$h_0(t)$
at the coalescence point is labelled. (b) Experimental data showing the dynamic shape of the interface in this plane. Notice that the darker-coloured markers correspond to earlier times and the lighter-coloured ones to later times. Markers are connected by a faint line that is intended to only serve as a guide for the eyes. (c) The interface profiles rescaled with
$h_0(t)$
. The black line is the self-similar profile in the
$xz$
plane. (d) Schematic of the interface in the
$yz$
plane, where
$a$
is the radius of a spherical cap. (e) Experimental data showing the dynamic shape of the interface in the
$yz$
plane. (f) The rescaled interface profiles with
$a=2.7$
mm.
Next, we present the interface profile at different times along the
$yz$
plane (figure 4
e). If the coalescence dynamics was similar to that of the Newtonian case in the thin-film regime, we would expect the shape of the interface along the
$yz$
plane to be parabolic. The parabolic profile is a limit of a circular segment shape when
$\theta \ll 1$
(Kaneelil et al. Reference Kaneelil, Pahlavan, Xue and Stone2022). Figure 4(f) shows the rescaled profiles in the
$yz$
plane and the black line is
$h(x=0,y,t)/h_0(t)=1- 1/2 (y/\sqrt {a h_0(t)} )^2$
. Note that the parameter
$a$
that appears in the rescaling is a geometric parameter that captures the outer length scale of the system (Kaneelil et al. Reference Kaneelil, Pahlavan, Xue and Stone2022). The data collapse well and show good agreement with the parabolic profile.
Thus, we have shown that the spatiotemporal dynamics of the interface in the
$xz$
and
$yz$
planes agrees with that expected for Newtonian coalescence, which suggests that the polymer has negligible effect on coalescence. We now show that the dynamic shape of the three-dimensional interface near the coalescence point can therefore be mapped to a two-dimensional self-similar curve (Kaneelil et al. Reference Kaneelil, Pahlavan, Xue and Stone2022). Figure 5(a) shows the time evolution of the three-dimensional interface for an experiment with 0.5 wt% PEO drops. Rescaling the height as
$S(\zeta ) = h(x,y,t)/[vt-y^2/(2a)]$
, where the similarity variable is
$\zeta =\theta x/[vt-y^2/(2a)]$
, collapses the data onto the self-similar shape, as shown in figure. 5(b).
Newtonian three-dimensional self-similarity also describes the coalescence of semi-dilute polymeric drops at small
$\theta$
. (a) Experimental data from the coalescence of 0.5 wt% PEO drops with
$\theta \approx 8.1^{\circ }$
, corresponding to
$De_\theta = 0.009$
and
$\textit{Ec}_\theta = 0.74$
, showing the three-dimensional shape of the interface near the coalescence point at early times (
$t= 0.05,\,0.15,\,0.22$
s). The darker-coloured markers correspond to earlier times and the lighter-coloured ones to later times. (b) Experimental data from the coalescence of 0.1, 0.5 and 1.0 wt% PEO drops at four different times and three different
$yz$
planes (total of 36 curves) rescaled according to the similarity solution. The rescaled data collapse onto the universal self-similar curve (black line).
In the small contact angle limit,
$\theta \ll 1$
, we showed experimentally that the coalescence dynamics of semi-dilute polymer solutions is unaffected by polymers, as predicted by our scaling analysis. However, we note that this is not a generic result for all thin-film problems. In the classical Landau–Levich–Derjaguin problem of drawing an object out of a bath of liquid and analysing the thickness of the liquid film left behind on the object, studies have shown that the thickness of the film is smaller for a weakly viscoelastic liquid relative to the Newtonian result (Datt, Kansal & Snoeijer Reference Datt, Kansal and Snoeijer2022), and larger when the viscoelastic effects are increased (Lee, Shaqfeh & Khomami Reference Lee, Shaqfeh and Khomami2002). The lack of an effect of the polymer that we observe can be attributed to the strong and unique dependence of the Deborah number on the contact angle,
$De_{\theta } \sim \theta ^3$
, as discussed in § 3.1. Notice that the scaling arises from the definition of the characteristic velocity in this problem
$u_c = \gamma \theta ^3/(3\mu )$
, which is a natural scale for capillary-driven thin-film flows. Thus, we would expect our result that thin-film flow dynamics is unaffected by polymers to hold for problems that are internally driven by capillarity and not externally driven, e.g. boundary driven, as in the case of the Landau–Levich–Derjaguin problem. This feature can be leveraged in applications involving thin-film flows of polymeric liquids, where it is advantageous to minimise elastic effects and attain outcomes similar to those achieved with Newtonian fluids.
4. The large contact angle limit
We perform two-dimensional numerical simulations of sessile drop coalescence using the open-source partial differential equation solver Basilisk C (Popinet Reference Popinet2009, Reference Popinet2015) using viscoelastic constitutive equations (Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018; López-Herrera et al. Reference López-Herrera, Popinet and Castrejón-Pita2019), along with some experiments to validate the numerical results in the large contact angle limit. We use the Oldroyd-B model to study the viscoelastic effects, since it is one of the simplest models, and show that it is sufficient to explain the regime of interest in this work.
While our experimental results show that coalescence is unaffected by polymer at small
$\theta$
, it has been shown experimentally to affect coalescence dynamics at larger
$\theta$
and when inertia is important (Varma et al. Reference Varma, Saha and Kumar2021; Dekker et al. Reference Dekker, Hack, Tewes, Datt, Bouillant and Snoeijer2022). Note that the characteristic scales that describe the flow dynamics are expected to be different here since inertia will not be negligible, whereas the
$\theta \ll 1$
approximation in the previous section naturally led to viscously dominated flows. The goal with the simulations is therefore to further probe the inertially dominated large-
$\theta$
limit to reveal the underlying mechanism that is leading to a departure from the Newtonian coalescence dynamics.
4.1. Dimensionless equations and numerical simulations
We rescale the governing equations from § 2 to make them dimensionless using the following characteristic scales that are relevant for inertial capillary-driven flows:
$u=\tilde {u}/U,\ t=\tilde {t}/\sqrt {\rho H^{3}/\gamma },\ y=\tilde {y}/H$
and
$x=\tilde {x}/H$
, where the typical velocity
$U$
is taken to be the capillary velocity
$U=H/t_c = (\gamma /\rho H)^{1/2}$
(Eggers Reference Eggers1997) and
$H$
is the height at the centre of the spherical cap-shaped drop. Notice that the characteristic scales are different from those in § 3.1, which was the viscously dominated thin-film regime. Non-dimensionalising (2.1b
) gives
\begin{equation} \frac {\textrm {D} \boldsymbol{u}}{\textrm {D} t} = -\boldsymbol{\nabla }p + \textit{Oh} {\nabla} ^2 \boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau^p}, \end{equation}
where the Ohnesorge number is
$\textit{Oh} = \mu_s /\sqrt {\rho \gamma H}$
, the pressure scale is chosen as
$p_c = \gamma /H$
, and the scale for the polymer stress was set to be the same as the pressure scale in order for polymer effects to appear in the first order. Non-dimensionalising (2.2) gives
\begin{gather} \boldsymbol{\tau }^p=Ec(\boldsymbol{A}-\boldsymbol{I}), \end{gather}
\begin{gather} \overset {\mathtt {\boldsymbol{\nabla }}}{\boldsymbol{A}}=-\frac {1}{De}(\boldsymbol{A}-\boldsymbol{I}), \end{gather}
\begin{gather} De\,\overset {\mathtt {\boldsymbol{\nabla }}}{\boldsymbol{\tau }^p} + \boldsymbol{\tau }^p = 2\, Ec De\, \boldsymbol{E}, \end{gather}
where the Deborah number and the elastocapillary number are defined as
\begin{equation} De = \frac {\lambda }{\sqrt {\rho H^{3}/\gamma }} \quad \textrm {and}\quad Ec=\frac {GH}{\gamma}. \end{equation}
Thus, the dimensionless parameters that affect the dynamics of the interface can be identified as
$\textit{Oh},\ De$
and
$\,Ec$
. In Basilisk C, the implementation of the Oldroyd-B model uses the so-called log conformation technique, where the logarithm of the conformation tensor
$\boldsymbol{A}$
is calculated (Fattal & Kupferman Reference Fattal and Kupferman2004; Hao & Pan Reference Hao and Pan2007; Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018; López-Herrera et al. Reference López-Herrera, Popinet and Castrejón-Pita2019). We used a modified version of the Oldroyd-B implementation that specifically identifies the modulus of the polymer instead of the solvent viscosity as one of the parameters, in addition to the relaxation time (Dixit et al. Reference Dixit, Oratis, Zinelis, Lohse and Sanjay2025).
The two-dimensional simulation is initialised in a square box of length
$L_0=2R$
with the shape of the interface defined as two symmetric circular segments, with all lengths rescaled by
$H$
, as follows:
\begin{align} h(x,t=0) &= {\mathcal H} (-x) \bigg [ \bigg ( \frac {1}{(1-\textrm {cos}(\theta ))^2} - (x+R)^2 \bigg )^{1/2} - \frac {1}{1-\textrm {cos}(\theta )} + 1\bigg ] \nonumber\\ &\quad+ {\mathcal H} (x) \bigg [ \bigg ( \frac {1}{(1-\textrm {cos}(\theta ))^2} - (x-R)^2 \bigg )^{1/2} - \frac {1}{1-\textrm {cos}(\theta )} + 1\bigg ] + h_{\infty }. \end{align}
Here,
$\mathcal H$
is the Heaviside step function,
$\theta$
is the angle that the circular segments make with the horizontal,
$R= [ ({1}/{(1-\textrm {cos}(\theta ))^2}) - ({\textrm {cos}(\theta )}/{(1-\textrm {cos}(\theta ))} )^2 ]^{1/2}$
is the rescaled radius of the base of the spherical cap and
$h_\infty$
is the arbitrarily small film thickness at the point where the two circular segments meet. Notice that
$h(x,t=0)$
is fully specified by a single parameter
$\theta$
. This definition of the interface shape eliminates the presence of a contact line. Figure 6(a) shows the simulation domain and the initial shape of the interface for a case with
$\theta =1.1\,\textrm {rad}$
and
$h_{\infty }=0.008$
(this value of
$h_{\infty }$
is kept constant for all simulations shown here and was verified to not affect the results). We impose a no-slip boundary condition on the bottom of the domain, symmetry conditions on the left and right and outflow on the top. The self-similar dynamics of the interface that we are interested in takes place at very early times of coalescence, when the height of the bridge is smaller than the macroscopic length scale of the drop
$h_0 \ll H$
, and should be unaffected by the symmetry conditions along the sides of the domain and other boundary effects.
Simulation set-up and comparison between simulation and experiments. (a) The initial shape of the interface
$h(x,t=0)$
, which follows (4.4), for a case with
$\theta =1.1$
and
$h_{\infty }=0.008$
is shown in the numerical domain of length
$L_0=2R$
. The inset shows the initial film thickness
$h_{\infty }$
at the coalescence point more clearly. (b) The temporal evolution of the height
$h_0(t)$
at the coalescence point. The filled markers are data from experiments and the open markers are data from simulations for
$\textit{Oh} = 0.004,\,\theta =1.1,\,De=0$
and
$Ec=0$
. (c) The experimental results (blue markers) overlaid on top of the simulation results (black lines) for
$\textit{Oh} = 0.004,\,\theta =1.1,\,De=0$
,
$Ec=0$
and
$t \approx [0,\,0.15]$
. (d) The shapes of the interface from experiment (dots) and simulation (lines) collapse onto a self-similar profile when rescaled with
$h_0(t)$
.
We use an adaptive quadtree grid where the maximum level of refinement
$L$
is specified, and the total number of grid points in the vertical or horizontal direction is
$2^L$
(Popinet Reference Popinet2015; van Hooft & Popinet Reference van Hooft and Popinet2022; Mostert, Popinet & Deike Reference Mostert, Popinet and Deike2022). We varied the grid refinement, as shown in Appendix A, to ensure that the results are grid-independent. We use
$L=11$
for all the simulations shown here, which corresponds to five grid points within a height of
$h_{\infty }$
.
4.2. Comparing experiments and simulations
We start by comparing the results between experiments and simulations for a Newtonian case with
$\textit{Oh} = 0.004,\,\theta =1.1,\,De=0$
and
$Ec=0$
, which corresponds to an experiment with water drops. Note that these large-
$\theta$
experiments are imaged from the side (see figure 2
a) using a high-speed camera (Phantom V2012; 29 000 fps). The large contact angles were achieved in experiments with glass microscope slides that were simply rinsed with water and dried with compressed air.
The results in figure 6(b) show the comparison of
$h_0/H$
as a function of time
$\tilde {t}/t_c$
between the simulations and experiments. There is excellent agreement between the rescaled simulated and experimental data. The slight deviation from the simulation results seen in the experimental data at very early times may be due to the error in measurement of the small
$h_0(t)$
values when it is close to the resolution of the camera.
In the Newtonian coalescence of sessile drops at large
$\theta$
and low
$\textit{Oh}$
, the dynamics of the height at the coalescence point follows
$h_0(t)/H = c_0 (\tilde {t}/t_c)^{2/3}$
. A value of
$c_0=0.89$
is used for the solid black line in figure 6(b), similar to previously reported results (Eddi, Winkels & Snoeijer Reference Eddi, Winkels and Snoeijer2013). Both the simulated and the experimental data in figure 6(b) are in agreement with this scaling.
Figure 6(c) shows the dynamic shape of the interface from the experiments with water (blue markers) overlaid on top of the results from the simulation (black lines) in rescaled units
$\tilde {x}/H$
and
$h/H$
. One artefact of the two-dimensional simulations is the more pronounced capillary waves seen at the interface leading to a wider interface shape than those observed in the inherently three-dimensional experiments (Keller, Milewski & Vanden-Broeck Reference Keller, Milewski and Vanden-Broeck2000; Eddi et al. Reference Eddi, Terwagne, Fort and Couder2008, Reference Eddi, Winkels and Snoeijer2013). Despite this slight discrepancy, the simulations and the experiments are in good agreement at early times of coalescence when there is a clear separation of scales. Figure 6(d) shows the shape of the interface obtained from both simulations (lines) and experiments (markers) rescaled according to the Newtonian self-similar scaling,
$h/h_0(t)$
and
$x/h_0(t)$
, which shows reasonable collapse.
We now use the numerical simulations to further study the dynamics of viscoelastic coalescence over a wide range of
$De$
and
$Ec$
in order to see how the coalescence dynamics deviates from Newtonian behaviour as a function of these dimensionless parameters.
4.3. The effect of
$De$
and
$Ec$
in the coalescence dynamics
The conspicuous effect of the polymer on the coalescence dynamics of large-
$\theta$
drops is seen in the shape of the interface, as shown in figure 1(b,c), where the curvature of the interface seems to be high in the presence of polymer. We now aim to understand the extent of this effect as a function of
$De$
and
$Ec$
. In all the following results, we keep the values of Ohnesorge number fixed at
$\textit{Oh}=0.01$
, ensuring that inertial–capillary effects dominate over viscous effects, the contact angle is
$\theta =1.1$
and focus solely on the physics as a function of
$De$
and
$Ec$
.
Figure 7(a) shows snapshots of the shape of the interface during coalescence at
$\tilde {t}/t_c=0.09$
for
$\textit{Oh}=0.01$
and
$\theta =1.1$
for various values of
$De$
and
$Ec$
. The coloured curve represents the shape of the interface at the given
$De$
and
$Ec$
, and the black curve represents the Newtonian case with
$De=Ec=0$
. When
$De$
and
$Ec$
are small, the shape of the interface overlaps with that of the Newtonian case. This limit corresponds to the case of dilute polymer solutions where the coalescence seems to be unaffected by polymers. As
$De$
and
$Ec$
increase, we see that the shape of the interface starts to change from that of the Newtonian response. The shape of the interface and therefore the viscoelastic effects seem to be more sensitive to
$Ec$
and show almost no change from the Newtonian results when
$Ec$
is sufficiently small.
Effect of
$De$
and
$Ec$
on the shape of the interface. (a) Snapshots of the interface at
$t = 0.09$
for
$\textit{Oh}=0.01$
and
$\theta =1.1$
for various
$De$
and
$Ec$
combinations. Note that the solid black lines correspond to the Newtonian limit with
$De=Ec=0$
, and the coloured lines correspond to the
$De$
and
$Ec$
shown on the axes. (b) The height at the coalescence point,
$h_0$
, is plotted as a function of time for various
$De$
and
$Ec$
. The solid black line corresponds to
$t^{2/3}$
while the dashed grey line marks
$t=0.09$
, the time at which the shapes in (a) are shown.
The results we see in the phase plane in figure 7(a) can be understood from the various limits of (4.2). First, for small
$De$
, i.e.
$De\rightarrow 0$
, we see from (4.2c
) that
$\boldsymbol{\tau }^p \approx 2 Ec\,De\,\boldsymbol{E}$
. If
$Ec\,De$
is also small, i.e.
$Ec \rightarrow 0$
, then we have
$\boldsymbol{\tau }^p \rightarrow 0$
, which is the Newtonian limit. If
$Ec\,De$
is finite or large, this scenario corresponds to the case where the total viscous response of the fluid is altered by the polymeric viscosity such that the deviatoric stress becomes
$(\textit{Oh}+Ec\,De){\nabla} ^2 \boldsymbol{u}$
. On the other hand, when
$Ec \rightarrow 0$
, we have
$\boldsymbol{\tau }^p \rightarrow \boldsymbol{0}$
independent of
$De$
, and the dynamics is purely Newtonian as seen in figure 7(a).
As
$De \rightarrow \infty$
, we have
$\overset {\mathtt {\boldsymbol{\nabla }}}{\boldsymbol{\tau }^p}=2Ec\boldsymbol{E}$
where the magnitude of
$Ec$
dictates whether the response is Newtonian or elastic. And finally as
$Ec \rightarrow \infty$
, we have
$\boldsymbol{\tau }^p \rightarrow \infty$
and
$\boldsymbol{E} \rightarrow \boldsymbol{0}$
corresponding to minimal deformation, as clearly seen by the undeformed interfaces in figure 7(a). We note that the actual results seen in experiments for
$De,\,Ec \rightarrow \infty$
might be quantitatively different from those seen in the simulations here since the Oldroyd-B model does not capture the finite extensibility of polymer chains.
The dynamics of the height of the bridge,
$h_0/H$
, is shown in figure 7(b) for several
$De{-}Ec$
combinations. Note that the data for
$Ec=10^3$
are not shown since the interface was effectively stationary. When
$De$
and
$Ec$
are sufficiently small, the temporal evolution follows the scaling
$t^{2/3}$
and the data are collapsed quite well with the Newtonian and inertial rescaling factors with the same prefactor of
$c_0=0.89$
(black line in figure 7
b). But, as the dimensionless parameters are increased, with
$De \approx 10^{-3}$
and
$Ec \approx 10^1$
and larger, a decrease from the
$t^{2/3}$
scaling is observed. Such a decrease in the scaling exponent as a function of the polymer concentration has been previously reported (Varma et al. Reference Varma, Saha and Kumar2021, Reference Varma, Dasgupta and Kumar2022a
; Rostami et al. Reference Rostami, Erb, Azizmalayeri, Steinmann, Stark and Auernhammer2025). However, we note that the scaling exponent in our study and in the results from Dekker et al. (Reference Dekker, Hack, Tewes, Datt, Bouillant and Snoeijer2022) goes lower than 0.5, perhaps since the way the drops come into contact is different from the droplet spreading method where the drops impact the substrate, spread, and contact each other (Varma et al. Reference Varma, Dasgupta and Kumar2022a
). We also point out that the departure from the Newtonian scaling can be seen in this study simply from the Oldroyd-B model without any need to consider shear-thinning and finite-extensibility behaviour. We now probe the evolution of the stress fields during coalescence to gain insight into this departure.
The larger curvature seen in the shape of the interface during viscoelastic coalescence (see Appendix C) suggests an increase in stress, presumably due to the contribution from
$\boldsymbol{\tau\!}^{p}$
, which we show is highly localised near the coalescence point. Figure 8(a) shows the dimensionless polymeric stress field (left-hand column) and the dimensionless inertial stress field (right-hand column) for
$De=10^{-2}$
and
$Ec=10^{-1}$
, while figure 8(b) shows the same for
$De=10$
and
$Ec=10$
. Note that we consider locally averaged stress fields, as explained in Appendix B, which are also independent of the grid size. In the former case, fluid inertia dominates and the effect of the polymeric stress is relatively minimal. In the latter case, however, polymeric stress dominates and is about an order of magnitude higher than the inertial stress. The region where this large polymeric stress occurs is also localised to near the middle because of the highly extensional flow created by the rising interface. Consequently, the larger stress must be balanced by a local increase in the pressure field, leading to a sharper interface profile.
The evolution of the polymeric and inertial stress fields for simulations with
$\textit{Oh}=0.01$
and
$\theta =1.1$
. The dimensionless polymeric stress field
$\tilde {\boldsymbol{\tau }}^p_{yy}/(\gamma /H)$
is plotted in the left-hand column and the dimensionless inertial stress
$\rho \tilde {u}_y^2/(\gamma /H)$
is plotted in the right-hand column, for four time steps. (a) A case with
$De=10^{-2},\,Ec=10^{-1}$
is observed to be dominated by inertia. (b) A case with
$De=10,\,Ec=10$
exhibits a much sharper interface shape and has a larger polymeric stress field. Note that the region of large polymeric stress is localised to the middle where the drops initially made contact.
Such a sharp shape of the interface was also observed at the rear end of gas bubbles rising in viscoelastic liquids, where the shape change was also attributed to locally large extensional stresses (Astarita & Apuzzo Reference Astarita and Apuzzo1965; Pilz & Brenn Reference Pilz and Brenn2007; Fraggedakis et al. Reference Fraggedakis, Pavlidis, Dimakopoulos and Tsamopoulos2016). These local extensional stresses caused by extended polymer molecules are generally present downstream of stagnation points and can be visualised in experiments using birefringence (Farrell & Keller Reference Farrell and Keller1978; Cressely, Hocquart & Scrivener Reference Cressely, Hocquart and Scrivener1979; Harlen, Rallison & Chilcott Reference Harlen, Rallison and Chilcott1990). Due to this local nature, large elastic stresses only occur close to these strands and the resulting flow can be modelled as a line distribution of forces in a Newtonian background fluid (Harlen et al. Reference Harlen, Rallison and Chilcott1990). In our study, the vertical rise of the interface during coalescence sets up a stagnation point at the substrate directly below the smallest height of the interface, namely
$h_0$
. At sufficiently high
$De$
and
$Ec$
, we expect the extensional flow resulting from coalescence to generate these strands of high elastic stresses, as seen in figure 8(b), and therefore alter the dynamics of the interface. We can probe the temporal evolution of the maximum values of the various stress fields to understand how the elastic stress develops as a function of our dimensionless parameters.
Figure 9 shows the evolution of the maximum value of the stress fields at early times of coalescence for various
$De$
and
$Ec$
. We keep
$\textit{Oh}=0.01$
and
$\theta =1.1$
as before. The blue marker represents max(
$\boldsymbol{\tau }^p_{yy}$
), red max(
$\boldsymbol{\tau }^{s}_{yy}$
), cyan max(
$\rho u_y^2$
) and green max(
$|p|$
). Note that the magnitude of pressure is taken since the pressure in the droplet phase will be negative during coalescence due to the reference pressure outside the drop being initially zero. Firstly, note that at very early times, the maximum value of pressure is independent of
$De$
and
$Ec$
since the curvature of the interface at the very initial moments of coalescence will likely only be a function of the contact angle
$\theta$
. Note also that the pressure here corresponds to the net stress in the system which is the same as the capillary stress since the interface shape reflects the state of stress. Once coalescence begins, the extent of how the other stresses balance pressure, or capillary stress, tells us about the dynamics of the interface. Quite remarkably, we see that at the very early moments of coalescence, the dominant balance is between capillary stress and inertia. Note that this will be true as long as
$\textit{Oh}$
is small and inertia dominates over viscous stress. When
$De$
and
$Ec$
are sufficiently small, inertial stresses primarily counterbalance the capillary stress throughout most of the coalescence process. However, as
$De$
and
$Ec$
increase, polymeric stress rapidly increases to become the predominant stress balancing capillary stress.
The maximum values of various stress fields as a function of
$\tilde {t}/t_c$
during the early times of coalescence for a wide range of
$De$
and
$Ec$
. The blue marker represents max(
$\boldsymbol{\tau }^p_{yy}$
), red max(
$\boldsymbol{\tau }^{s}_{yy}$
), cyan max(
$\rho u_y^2$
) and green max(
$|p|$
). Note that all stresses are rescaled by
$\gamma /H$
. The dashed black line represents
$h_0^{-1} \propto t^{-2/3}$
which is the rate at which capillary and inertial stress decays in the absence of polymers.
In the initial stages of coalescence before the polymeric stress takes over, the capillary stress and the balancing inertial stress are reasonably expected to decay independent of
$De$
and
$Ec$
as
$h_0^{-1} \propto t^{-2/3}$
(black dashed line in figure 9). As polymeric stress increases to become the dominant contribution, the rate of decay of the capillary stress decreases, shown by the departure from the
$t^{-2/3}$
behaviour in figure 8 (see green data points). This decrease in the rate of decay of the capillary stress directly leads to a slower growth of
$h_0$
. Furthermore, since the curvature of the interface is tied to the stress field, a greater curvature is observed in the shape of the interface (see figures 7 and 8). This increase in the polymeric stress, which is spatially local in nature, and the subsequent slowing in the rate of decay of the capillary stress are what cause the coalescence dynamics to deviate from the inertial and Newtonian behaviour.
Interestingly, near
$Ec \approx 10$
when the departure from the Newtonian dynamics is observed, the polymeric stress still requires a finite amount of time to grow and become the dominant contribution, allowing for the capillary and the inertial stresses to set the early-time dynamical behaviour. This finite time scale of the growth of the polymer stress is what allows the Newtonian scaling of the bridge height
$h_0(t) \sim t^{2/3}$
to persist. However, we see that the dynamics slowly deviates from the
$t^{2/3}$
scaling over time (see figure 7
b) once the polymeric stress has become the dominant contribution. While the departure from the Newtonian dynamics of the temporal evolution of the bridge is subtle and gradual, the effect on the shape of the interface is effectively instantaneous, since the shape of the interface must directly reflect the instantaneous state of stress. It is also more noticeable because the polymeric stress is highly localised (see figure 8). This is clearly observed when analysing the temporal evolution of the curvature of the interface (see Appendix C), which shows that the curvature for the viscoelastic case is much larger than its Newtonian counterpart whereas the temporal scaling remains the same.
Increasing
$Ec$
seems to decrease the time scale of the growth of the polymer stress, leading to an earlier transition to polymeric stress dominance and allowing the polymeric stress to reach a greater magnitude. As a result, the interface shape is altered more significantly with larger
$Ec$
. Increasing
$De$
increases the time scale over which the polymeric stress starts to decay and essentially dictates whether or not the transition to polymeric stress dominance occurs during the early stages of coalescence. Thus, we show that elasticity affects coalescence dynamics only at sufficiently large
$De$
and
$Ec$
due to the extensional flow giving rise to locally large elastic stress strands. We believe that more experimental and theoretical studies are necessary to fully understand this complicated process of viscoelastic coalescence and reveal the self-similar dynamics which must account for the transition we see as a function of
$De$
and
$Ec$
. Future work may also probe the effect of varying
$\textit{Oh}$
and
$\theta$
and consider more sophisticated viscoelastic constitutive laws.
5. Conclusions
The coalescence of polymer drops was explored using experiments in the context of small
$\theta$
and using experiments and direct numerical simulations in the context of large
$\theta$
. For the small-
$\theta$
case, we imaged the three-dimensional shape of the interface during coalescence using FS-SS and found that the coalescence dynamics is the same as that of viscous Newtonian drops. The negligible effect of the polymer on coalescence was attributed to the low Deborah number of the system, which scales with
$\theta$
as
$De \sim \theta ^3$
in the small-
$\theta$
limit. When
$\theta \ll 1$
, we showed that the effect of the polymers only appears at
${\mathcal O}(\theta )$
and the coalescence dynamics simply follows the Newtonian self-similar dynamics for the thin-film regime. This result is important because viscoelasticity is known to affect the thin-film dynamics in other problems, such as the Landau–Levich–Derjaguin problem. Yet, our results show that coalescence is a scenario where the elastic effects can be effectively disregarded within the thin-film flow regime.
In the large-
$\theta$
case where inertia dominates the dynamics, coalescence is altered by the presence of polymers. We used direct numerical simulations with the Oldroyd-B model to study this limit. The two parameters that capture the effect of viscoelasticity were the Deborah number
$De$
and the elastocapillary number
$Ec$
. As
$De$
and
$Ec$
are increased, we showed that the polymeric stress becomes the dominant stress contribution balancing capillary pressure. Revealing the stress fields inside the drop and the temporal evolution of the various stress fields, we provide more insight into how the coalescence dynamics starts to deviate from the classical Newtonian dynamics for viscoelastic drops.
Acknowledgements
We thank D.B. Shaw and J. Wu for helpful discussions regarding the numerical simulations. We thank V. Sanjay for discussions regarding a modified Oldroyd-B implementation for Basilisk C.
Funding
This work was supported by the United States National Science Foundation grant CBET-2242512 to L.D. and CBET-2246791 to H.A.S.
Data availability statement
The code used for the simulations in this study are available at https://github.com/paulkaneelil/viscoelasticCoalescence.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Grid refinement
We varied the grid refinement to ensure that the results are grid-independent. Figure 10 shows the result from a grid-independence study in which
$\textit{Oh} = 0.09,\,\theta =1.1,\,De=3.6$
and
$Ec=0.04$
. Figure 10(a) shows the temporal evolution of the height
$h_0(t)$
at the coalescence point for three different
$L$
values. Note that the lengths are rescaled by
$H$
and the time is rescaled by the inertial–capillary time:
$t_c = \sqrt {\rho H^{3}/\gamma }$
. Figure 10(b) shows the shape of the interface at different time steps for different
$L$
values. The overlapping shape of the interface in figure 10 suggests that the data are well resolved at this grid refinement.
Results from varying the grid size in the simulations with
$\textit{Oh} = 0.09,\,\theta =1.1,\,De=0.01$
and
$Ec=0.04$
. (a) The temporal evolution of the height
$h_0(t)$
at the coalescence point is shown for three different grid refinements:
$2^{9},\,2^{10}$
and
$2^{11}$
. The solid black like shows the
$t^{2/3}$
scaling. (b) The shape of the interface for the three different grid refinements at three different time points.
Appendix B. Calculating average stress
All the stress fields from the simulations were locally averaged over a square of size
$\epsilon =0.007$
to make the results independent of the grid size. The average stress can be represented as
\begin{equation} \sigma = \frac {1}{\epsilon ^2} \iint \hat {\sigma }\,\mathrm{d}S, \end{equation}
where
$\mathrm{d}S$
is a square of length
$\epsilon$
. The value of
$\epsilon$
was varied and was taken to be such that the result is independent of
$\epsilon$
itself, and turned out to be slightly larger than the length of a
$3 \times 3$
grid.
Appendix C. Interfacial curvature
The interfacial curvature in the
$xz$
plane near the coalescence point (
$x=0$
) was calculated to quantify the sharpness of the interface. Given the shape of the interface
$h(x,t)$
, the curvature is defined as
\begin{equation} \kappa = \frac {h_{xx}}{(1+h_x^2)^{3/2}} \approx h_{xx}, \end{equation}
where the approximation is made since the slope vanishes at the coalescence point for symmetric coalescence. The curvature was calculated by fitting a parabola to the interface data,
$h(x,t) = h_0(t) + ({1}/{2})(\kappa (t)\,x^2)$
, considering the appropriate amount of data points that minimised the root-mean-square error of the fit. The results for the experiments are shown in figure 11(a) and for the simulations are shown in figure 11(b). Both the experimental and simulated data clearly show an increase in interfacial curvature in the presence of polymers. This quantification confirms the observation that was made in figure 1(b,c).
The interfacial curvature for Newtonian drops scales as
$\kappa \sim 1/h_0 \sim t^{-2/3}$
, which is the black line in figure 11. Note that as
$De$
and
$Ec$
increase, in figure 11(b), the data deviate away from this scaling.
The time evolution of the interfacial curvature near the coalescence point. (a) Curvature calculated from the experimental results where the blue circles correspond to a Newtonian case with
$\textit{Oh} = 0.004,\theta =1.1,\,De=0$
and
$Ec=0$
(experiment with water drops) and the red circles correspond to a non-Newtonian case with
$\textit{Oh} = 0.09,\,\theta =1.1,\,De=3.6$
and
$Ec=0.04$
(experiment with 0.5 wt% PEO drops). The black line represents
$t^{-2/3}$
. (b) Curvature calculated from the simulation results. The data have the same legend as that shown in figure 7 in the main text. The black line represents
$t^{-2/3}$
.
Appendix D. Free-surface synthetic schlieren: image processing
The FS-SS utilises the fact that the curved shape of the drop acts as a lens to distort the image of the dot pattern that is placed under the drop. By measuring the distortion, or the displacement field of the dot pattern compared with a reference image, we can calculate the shape of the interface that caused the distortion. The displacement field is calculated using PIVLab (Thielicke & Stamhuis Reference Thielicke and Stamhuis2014). Figure 12(a) shows a reference image of the dot pattern above which the coalescence experiment is performed. The dot pattern is produced by evaporating a dense suspension of polydispersed polystyrene particles. A microscope slide of thickness 1 mm is placed above the dot pattern, with approximately 0.3 mm space in between. Thus, the distance between the dot pattern and the plane of coalescence is
$\Delta h = 1.3$
mm. The value of
$\Delta h$
is used in image processing to calculate the interface profile. Figure 12(b) shows the distorted image of the dot pattern at
$t=1$
s into the coalescence of 1 wt% PEO drops, and figure 12(c) shows the displacement field vectors for this frame. Note that not all the calculated vectors are shown in the image to avoid crowding.
Image processing to derive height data from FS-SS imaging. (a) Image showing the reference configuration of the dot pattern without a drop on top. (b) Snapshot of coalescence of 1 wt% PEO drops at
$t=1$
s (
$\theta \approx 11.7^{\circ }$
). Notice that the image is slightly distorted compared with the reference image. (c) Displacement field from PIVLab for the
$t=1$
s image showing the extent of the distortion caused by the drop on top of the dot pattern. All scale bars represent 0.5 mm. (d) Unfiltered and (e) filtered data showing the height evolution at the coalescence point. (e) Reconstructed three-dimensional profile of the interface for
$t=1$
s.
The displacement field is converted to the height profile of the interface using functions in the PIVMat toolbox in MATLAB (Moisy Reference Moisy2025), which is based on the surface height reconstruction method delineated in Moisy et al. (Reference Moisy, Rabaud and Salsac2009). Figure 12(d) shows the height
$h_0(t)$
at the coalescence point from the reconstructed interface profile data. The spatial resolution of our FS-SS experiments is of the order of
$0.1$
μm and is a function of the imaging window size and the resolution of the camera. However, note that the data show high-frequency oscillations. This noise is due to the internal fan of the high-speed camera which produces vibrations that are large enough to be observed within the spatiotemporal resolution with which we operate. We use a low-pass filter to attenuate the high-frequency oscillations, and the filtered data are shown in figure 12(e). No further smoothing is done on the experimental data. Figure 12(f) shows the reconstructed three-dimensional shape of the interface corresponding to
$t=1$
s.
Open access
