# Transition to elasto-capillary thinning dynamics in viscoelastic jets

Konstantinos Zinelis 1,2  Thomas Abadie 1,3 Gareth H. McKinley 2  Omar K.
Matar1<EMAIL_ADDRESS>1Department of Chemical Engineering, Imperial
College London, London SW7 2AZ, UK 2Department of Mechanical Engineering,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA 3School of
Chemical Engineering, University of Birmingham, Birmingham B15 2TT, UK

###### Abstract

We perform simulations of an impulsively-started, axisymmetric viscoelastic
jet exiting a nozzle and entering a stagnant gas phase using the open-source
code Basilisk. This code allows for efficient computations through an
adaptively-refined volume-of-fluid technique that can accurately capture the
deformation of the liquid-gas interface. We use the FENE-P constitutive
equation to describe the viscoelasticity of the liquid and employ the log-
conformation transformation, which provides stable solutions for the evolution
of the conformation tensor as the jet thins down under the action of
interfacial tension. For the first time, the entire jetting and breakup
process of a viscoelastic fluid is simulated, including the pre-shearing flow
through the nozzle, which results in an inhomogeneous initial radial stress
distribution in the fluid thread that affects the subsequent breakup dynamics.
The evolution of the velocity field and the elastic stresses in the nozzle are
validated against analytical solutions where possible, and the early-stage
dynamics of the jet evolution are compared favourably to the predictions of
linear stability theory. We study the effect of the flow inside the nozzle on
the thinning dynamics of the viscoelastic jet (which develops distinctive
“beads-on-a-string” structures) and on the spatio-temporal evolution of the
polymeric stresses in order to systematically explore the dependence of the
filament thinning and breakup characteristics on the initial axial momentum of
the jet and the extensibility of the dissolved polymer chains.

##  1 Introduction

Spray formation, which involves the disintegration of a continuous liquid
stream as it enters into a stagnant gaseous phase, is an important aspect of
many industrial and biological processes (Villermaux, 2007). Some
representative examples include inkjet printing processes (Basaran et al.,
2013; Lohse, 2022), the dispersal of fertilizers and pesticides on plants (Xu
et al., 2021), as well as human sneezing (Scharfman et al., 2016). In
particular, sprays are the result of an atomisation process in which a liquid
jet is destabilized and undergoes breakup into ligaments, threads, and
eventually droplets, under the action of capillary, inertial, and viscous
forces. The breakup process involves complex interfacial topological
transitions featuring pinch-off singularities where the filament radius
locally goes to zero.

The addition of polymers to a Newtonian solvent endows the resultant polymeric
solution with elasticity, which can significantly influence these
destabilization and breakup processes; due to the increased extensional
resistance to elongation that arises as a result of stretching of the
polymeric chains. Examples of viscoelastic systems include paints, inks,
industrial thickeners, anti-misting polymer agents, and human saliva or mucus.
Achieving a fundamental understanding of the fragmentation process in the
presence of viscoelastic effects will facilitate optimisation and control of
the droplet size distribution associated with sprays of polymeric solutions.

In this paper, we first consider the phenomenology of the thinning and breakup
of thin filaments of a polymeric solution that leads to beads-on-a-string
(BOAS) structures, for which there is no analogue in simple fluids; these
structures correspond to a series of almost cylindrical filaments connecting
spherical beads (Clasen et al., 2006). Extensive experimental and numerical
work has been conducted to understand the processes leading to the formation
of BOAS structures in a thinning viscoelastic filament. The typical
formulation makes use of slender jet profile approximation, resulting in a
one-dimensional description (Clasen et al., 2006; Eggers & Villermaux, 2008)
of the jet radius, axial velocity, and polymeric stress components. More
recently, the self-similar profiles of a viscoelastic thread during the
thinning process have also been computed with Direct Numerical Simulations
(DNS) (Turkoz et al., 2018; Snoeijer et al., 2020), highlighting the local
importance of the polymeric shear stress components to the final pinch-off.

A characteristic exponential rate-of-thinning in the viscoelastic filament, in
which the initially capillary-driven deformation of the fluid interface is
eventually balanced by fluid elasticity, has also been observed both
experimentally (Bazilevesky et al., 1990; Entov & Hinch, 1997; Amarouchene et
al., 2001; Clasen et al., 2006; Deblais et al., 2020) and through numerical
simulations (Bousfield et al., 1986; Étienne et al., 2006; Bhat et al., 2008,
2010; Morrison & Harlen, 2010; Turkoz et al., 2018; Eggers et al., 2020). The
analytic studies have primarily considered infinitely-extensible polymeric
chains that can be described by the Oldroyd-B (Bird et al., 1987) constitutive
equation (Li & Fontelos, 2003; Clasen et al., 2006; Ardekani et al., 2010),
while computational simulations have investigated the final breakup of the
thread that results when finite extensibility of macromolecules is
incorporated by using the FENE-P (Bird et al., 1987) model to describe the
polymeric stress evolution (Anna et al., 2001; Fontelos & Li, 2004; Wagner et
al., 2005). In particular, the role of finite extensibility in controlling the
final breakup time (or length) of viscoelastic threads has been examined, and
an analytical solution, which describes the local pinch-off dynamics, has been
derived (Entov & Hinch, 1997; McKinley, 2005; Wagner et al., 2015). In
addition, the generation of satellite droplets, as well as the influence of
viscosity and the fluid relaxation time on a viscoelastic filament initially
at rest have been studied via two-dimensional simulations (Liu et al., 2023).
Recently, numerical simulations have also been performed of viscoelastic jets
and single droplet breakup in an inkjet-printing configuration (without
considering the flow in a nozzle) using the Oldroyd-B model (Turkoz et al.,
2021), and an adaptive remeshing algorithm was introduced to resolve the time-
dependent evolution of very thin viscoelastic threads.

To understand the dynamical behaviour and the consequences of the viscoelastic
nature of non-Newtonian liquids, we must quantify the extensional rheological
properties of this class of materials. To this end, various experimental
protocols accounting for different initial liquid configurations have been
developed. The Capillary Breakup Extensional Rheometer (CaBER) (Bazilevesky et
al., 1990; Yesilata et al., 2006) for various dilute and semi-dilute polymeric
solutions, and the Rayleigh Ohnesorge Jetting Extensional Rheometer (ROJER)
(Keshavarz et al., 2015) are the most frequently applied methods to measure
the extensional rheology of viscoelastic fluids. The latter technique is
equivalent to a “flying” CaBER, and allows the examination of a wider range of
fluid relaxation times. Additional results for both industrial and biological
fluids have also been reported using other recently developed instruments that
exploit capillary-driven breakup such as Dripping-onto-Substrate (DoS)
Rheometry (Dinic et al., 2015; Dinic & Sharma, 2019; Martínez Narváez et al.,
2021; Lauser et al., 2021) as well as for more standard dripping experiments
which focus on the transition to elasticity-dominated thinning (Amarouchene et
al., 2001; Wagner et al., 2005; Rajesh et al., 2022).

Although the establishment of an elasto-capillary (EC) balance that results in
an exponential rate-of-thinning in the fluid thread has been validated in each
of these extensional rheometry configurations, recent experimental and
numerical considerations have argued that there can be systematic differences
in the local rate of thinning observed in CaBER and ROJER experiments (Mathues
et al., 2018). This may be due to non-zero initial values of the axial
stresses in the filament which develop as the liquid jet is expelled through
the nozzle exit in the ROJER configuration, subtly altering the tensile force
balance on a thin viscoelastic filament and leading to a $33\%$ faster
exponential decrease of the local jet radius. These faster dynamics were
reported for the first time in the limit of very low jet ejection flowrates
due to the so-called “gobbling phenomenon”, with the conventional elasto-
capillary balance expected in the CaBER instrument being established for
larger flowrate values (Keshavarz et al., 2015; Sharma et al., 2015). The
origins of these different thinning dynamics and their dependence on the
viscoelastic properties of the fluid and process parameters such as the nozzle
radius and the jet velocity remain poorly understood.

The present work aims to address the issues highlighted above and facilitate
the design of robust extensional rheometric instrumentation for accurately
measuring the rheological properties of weakly elastic complex fluids, such as
the relaxation time, the transient extensional viscosity, and the
corresponding strain rate. This requires developing a quantitative
understanding of the delicate interplay among viscous, inertial, capillary,
and elastic forces in low-speed axisymmetric jets in order to establish robust
foundations for future studies of the complex dynamics of viscoelastic sprays
and the ensuing droplet size distributions that develop at higher flowrates.
To achieve this, axisymmetric numerical simulations are carried out over a
wide range of system parameters using the open-source code Basilisk (Popinet,
2009) which incorporate the viscoelastic shearing flow upstream of the nozzle
exit in addition to the subsequent capillarity-driven evolution of the jet and
the formation of beads-on-a-string morphologies.

The rest of this article is organised as follows. In Section 2, the problem
formulation and numerical procedure used to carry out the computations are
outlined. A discussion of the results is provided in Section 3 highlighting
the stress profiles that develop due to the flow within the nozzle and the
subsequent evolution of the jet towards breakup following its exit from the
nozzle. Careful attention is paid to the evolution of the jet thinning
characteristics with changes in the flowrate and the finite extensibility of
the dissolved polymer. Finally, concluding remarks are provided in Section 4.

##  2 Formulation and Methodology

We first provide details of the problem formulation which encompasses the flow
configurations, the governing equations including the constitutive relation
used to describe the fluid viscoelasticity, and the numerical methodology
deployed to carry out the simulations.

###  2.1 Governing equations and numerical method

The simulation setup for an axisymmetric jet of an incompressible fluid of
density $\rho_{l}$ issuing from a nozzle of length $\ell_{nozzle}$ and initial
radius $R_{0}$, is presented in Figure 1. The fluid corresponds to a
viscoelastic polymer solution whose polymeric chains have finite extensibility
$L$; the total fluid zero-shear viscosity is $\eta_{0}=\eta_{p}+\eta_{s}$
wherein $\eta_{p}$ and $\eta_{s}$ denote the polymer and solvent contributions
to the viscosity, respectively. Additionally, $\beta$ is the fluid viscosity
ratio $\eta_{s}/\eta_{0}$ which determines the relative polymeric and solvent
contributions to the total dynamic viscosity of the polymer solution. The jet
is surrounded by a gas of density $\rho_{g}$ and viscosity $\eta_{g}$, and
gravitational effects are neglected in the present work. In the simulations,
we take the characteristic liquid viscosity to be $\eta_{l}\equiv\eta_{0}$,
where $\eta_{l}$ stands for the viscosity of the liquid phase. A pressure
gradient (see section 2.2 below) is imposed along the nozzle in the axial
direction $x$ (here, the considered axisymmetric coordinates are
$(r,\theta,x$)) leading to the development of a parabolic velocity profile
with mean velocity $U_{0}$ within the nozzle. As shown schematically in Figure
1, following its exit from the nozzle, the jet evolves downstream with a
characteristic wavelength $\lambda_{w}$ due to the imposed perturbations,
undergoing an increasingly pronounced deformation over a length $\ell_{jet}$
until ultimately a breakup event occurs. Jet breakup proceeds via the
development of a thin thread that connects the leading drop with the rest of
the liquid core and polymer stresses in this highly stretched thread act to
retard its eventual detachment. Following a fluid element of fixed Lagrangian
identity in this thinning ligament reveals an exponential decrease in the
radius with time. As will be discussed below, this exponential local decrease
of $R_{min}(t)$ results from the establishment of a local elasto-capillary
balance in the thinning thread.

Figure 1: Schematic representation of the flow depicting the jetting and
eventual breakup of a viscoelastic fluid issuing from a nozzle surrounded by a
gaseous phase. The flow upstream of the nozzle exit is fully-developed, driven
by an applied pressure gradient. The solution domain is highlighted in red.

The jet dynamics are governed by the one-fluid formulation of the continuity
and momentum equations, which are respectively expressed by:

$\boldsymbol{\nabla\cdot\mathbf{u}}=0,$ (1)
$\rho\left(\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\boldsymbol{\nabla}\mathbf{u}\right)=-\boldsymbol{\nabla}p+\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}+\gamma\kappa\mathbf{n}\delta_{s},$
(2)

where $t$, $\rho$, $\mathbf{u}$, $p$, $\boldsymbol{\sigma}$, $\gamma$,
$\kappa$, $\mathbf{n}$, and $\delta_{s}$ stand for time, local density,
velocity, pressure, the total stress, (constant) surface tension, interfacial
curvature, the outward-pointing unit vector to the interface, and the Dirac
delta function (zero everywhere except at the interface), respectively. Here,
$\gamma\kappa{\bf n}\delta_{s}$ denotes the surface tension forces distributed
in the cells in the vicinity of the interface with the Continuum Surface Force
method (Popinet, 2009, 2018). Given the viscoelastic nature of the fluid, the
total stress is defined as the sum of the solvent and polymeric stresses,
$\boldsymbol{\sigma}=\boldsymbol{\sigma}_{s}+\boldsymbol{\sigma}_{p}$, where
$\boldsymbol{\sigma}_{s}=\eta_{s}(\boldsymbol{\nabla}\mathbf{u}+(\boldsymbol{\nabla}\mathbf{u})^{T})$
is the viscous contribution to the total stress tensor, and
$\boldsymbol{\sigma}_{p}$ is the polymeric stress tensor defined in the
present work by the FENE-P constitutive equation:

$\boldsymbol{\sigma}_{p}=\frac{\eta_{p}}{\tau}\left(\frac{\mathbf{A}}{1-\frac{\operatorname{tr}(\mathbf{A})}{L^{2}}}-\mathbf{I}\right);$
(3)

here $\tau$ is the single characteristic relaxation time of the viscoelastic
fluid, $L^{2}$ provides a measure of the finite extensibility of the polymeric
chains, and $\mathbf{A}$ is the dimensionless conformation tensor of the
finitely extensible nonlinearly elastic (FENE) dumbbells that model the
evolution of the polymer configuration, which is governed by the following
equation:

$\frac{\partial\mathbf{A}}{\partial
t}+\mathbf{u}\cdot\boldsymbol{\nabla}\mathbf{A}-\left(\boldsymbol{\nabla}\mathbf{u}\cdot\mathbf{A}+\mathbf{A}\cdot\boldsymbol{\nabla}\mathbf{u}^{T}\right)=-\frac{1}{\tau}\left(\frac{\mathbf{A}}{1-\frac{\operatorname{tr}(\mathbf{A})}{L^{2}}}-\mathbf{I}\right).$
(4)

The above equations are rendered dimensionless by using the nozzle radius
$R_{0}$, the Rayleigh velocity
$U_{R}=\sqrt{\gamma/\left(\rho_{l}R_{0}\right)}$, the Rayleigh time scale
$t_{R}=R_{0}/U_{R}=\sqrt{\rho_{l}R_{0}^{3}/\gamma}$, and the capillary
pressure $\rho_{l}U_{R}^{2}=\gamma/R_{0}$ as the characteristic length,
velocity, time, and pressure/stress scales, respectively. Introduction of this
scaling into Eqs. (1)-(4) leads to the following dimensionless equations:

$\boldsymbol{\tilde{\boldsymbol{\nabla}}\cdot\mathbf{\tilde{u}}}=0,$ (5)
$\tilde{\rho}\left(\dfrac{\partial\mathbf{\tilde{u}}}{\partial\tilde{t}}+\mathbf{\tilde{u}}\cdot\tilde{\boldsymbol{\nabla}}\mathbf{\tilde{u}}\right)=-\tilde{\boldsymbol{\nabla}}\tilde{p}+Oh\left(\beta\tilde{\boldsymbol{\nabla}}\cdot\boldsymbol{\tilde{\sigma}}_{s}+\frac{\left(1-\beta\right)}{De}\tilde{\boldsymbol{\nabla}}\cdot\boldsymbol{\tilde{\sigma}}_{p}\right)+\tilde{\kappa}\tilde{\delta}\mathbf{n},$
(6)
$\frac{\partial\mathbf{A}}{\partial\tilde{t}}+\mathbf{\tilde{u}}\cdot\boldsymbol{\tilde{\boldsymbol{\nabla}}}\mathbf{A}-\left(\boldsymbol{\tilde{\boldsymbol{\nabla}}}\mathbf{\tilde{u}}\cdot\mathbf{A}+\mathbf{A}\cdot\boldsymbol{\tilde{\boldsymbol{\nabla}}}\mathbf{\tilde{u}}^{T}\right)=-\frac{1}{De}\left(\frac{\mathbf{A}}{1-\frac{\operatorname{tr}(\mathbf{A})}{L^{2}}}-\mathbf{I}\right),$
(7)

where $De=\tau/(R_{0}/U_{R})=\tau\sqrt{\gamma/\rho_{l}R_{0}^{3}}$ denotes the
Deborah number, which represents the ratio of the polymer relaxation time to
the Rayleigh time scale, and $Oh=\eta_{l}/\sqrt{\rho_{l}\gamma R_{0}}$ is the
Ohnesorge number that reflects the competition between capillary, inertial,
and viscous forces. The tildes designate dimensionless variables.

Additionally, the importance of the polymer elasticity in the neck can also be
understood through a local strain rate
$\dot{\epsilon}_{min}=-2D\left(\log(R_{min}^{[\alpha]})\right)/Dt$ for each
local minimum radius observed in each neck that is formed and subsequently
develops into a thin cylindrical ligament between two consecutive primary
beads (as depicted in Figure 1), where the primary beads are labelled as
$[\alpha]=$ A,B,C,… (starting from the beads furthest from the nozzle). When
scaled with $\tau$, this local strain rate in a material element, as it is
advected downstream, corresponds to the local Weissenberg number in the
thinning ligament:

$Wi=-2\tau\frac{D\log(R_{min}^{[\alpha]})}{Dt}.$ (8)

This definition will be used to characterise the local rate of jet thinning in
each neck that develops along the corrugated jet. When capillarity and
elasticity locally govern the dynamics in the elasto-capillary regime as
breakup is approached it is expected that the Weissenberg number will approach
a constant value $Wi\rightarrow 2/3$ (Entov & Yarin, 1984; Bazilevesky et al.,
1990; Entov & Hinch, 1997; Tirtaatmadja et al., 2006).

As a result of the coil-stretch transition, the local polymeric stresses in a
fluid element typically exhibit a steep increase under the influence of
extensional deformations at $Wi\geq 0.5$. In such configurations, numerical
challenges can emerge during the computation of the polymeric stress tensor
and its divergence which is required for Eq. (3). This is known as the High-
Weissenberg Number Problem (HWNP) (Renardy, 2000). To circumvent the
occurrence of any undesired numerical instability, the log-conformation
transformation (Fattal & Kupferman, 2005) is employed which introduces the
logarithmic function $\boldsymbol{\Uppsi}$ of the conformation tensor
$\mathbf{A}$, which can be computed via a diagonal transformation such that:

$\mathbf{A}=\mathbf{R}\ \boldsymbol{\Uplambda}\ \mathbf{R}^{T},$ (9)
$\boldsymbol{\Uppsi}=\log\mathbf{A}=\mathbf{R}\ \log\boldsymbol{\Uplambda}\
\mathbf{R}^{T}$ (10)

where $\boldsymbol{\Uplambda}$ is the diagonal matrix containing the principal
eigenvalues of $\mathbf{A}$.

A Volume-of-Fluid (VOF) approach (Popinet, 2009) is used to capture the
deforming jet interface by advecting the volume fraction $c$ of the liquid
phase in every computational cell; the advection equation for $c$ is given by

$\frac{\partial
c}{\partial\tilde{t}}+\mathbf{\tilde{u}}\cdot\boldsymbol{\tilde{\boldsymbol{\nabla}}}c=0.$
(11)

Using the VOF formulation, which corresponds to a one-fluid approach to solve
the two-phase flow, the density $\rho$ and viscosity $\eta$ are then
respectively given by

$\displaystyle\tilde{\rho}$ $\displaystyle=c+(1-c)\dfrac{\rho_{g}}{\rho_{l}},$
(12) $\displaystyle\tilde{\eta}$
$\displaystyle=c+(1-c)\dfrac{\eta_{g}}{\eta_{l}};$ (13)

here, the characteristic density is chosen to be $\rho_{l}$, and the tilde
decoration is suppressed henceforth for brevity. The open-source code
_Basilisk_ (Popinet, 2009; Turkoz et al., 2018; López-Herrera et al., 2019;
Turkoz et al., 2021) is used here to carry out the computations. A piecewise
linear interface calculation (PLIC) technique is used for reconstructing the
interface (Popinet, 2009; López-Herrera et al., 2019). The surface tension-
dominated flow disintegration of a liquid jet is modelled with high accuracy
thanks to the well-balanced numerical discretisation combined with the height-
function method to calculate the geometrical properties of the interface
(Popinet, 2009, 2018).

###  2.2 Numerical set-up

The simulation domain for the axisymmetric simulations of the jet is a plane
of dimensions $100R_{0}\times 100R_{0}$ as shown in Fig. 1. The bottom
boundary is the axial symmetry axis, while zero-gradient Neumann boundary
conditions are imposed at the left and right boundaries for all the velocity
and polymeric stress components. A pressure gradient is imposed at the left
boundary to drive the fluid through the nozzle while atmospheric pressure is
imposed at the right boundary through a Dirichlet condition. For the initial
conditions of the jetting process at $t=0$, we consider the liquid initially
at rest in the region $0\leq r\leq R_{0}$ and $0\leq x\leq\ell_{nozzle}$, and
the dissolved polymers in the fluid to be initially unstretched in the inlet
of the nozzle, with $A_{xx}=A_{rr}=1$.

An adaptive mesh refinement (AMR) technique (Popinet, 2003) following the
quadtree-like structure available in _Basilisk_ is used to refine the cells
based on the location of the interface and the nozzle, as well as in the
regions where large gradients of the axial component of the polymeric stress
occur. Specifically, starting from a base grid resolution of $8\times 8$
square cells for the entire domain and refining up to an initial minimum cell
size of $\Delta x_{minimum}=0.09$ around the nozzle region, the adaptive
scheme refines up to three maximum levels of refinement, LVL, gradually
increasing from LVL=12 to LVL=14, which corresponds to a minimum square cell
size of $\Delta x_{minimum}=0.02$ and $\Delta x_{minimum}=0.006$,
respectively; this provides sufficient resolution to simulate the dynamics
accurately as the pinch-off of the filament is approached. More details of the
mesh convergence study are provided in Appendix A.

To simulate the flow within the nozzle, a simplified variant of the Immersed
Boundary Method (IBM) (Aniszewski et al., 2020) is employed to model the
velocity field with no-slip and no-penetration conditions imposed at the solid
walls. The parameter values used in the simulations are provided in Table 1.
Starting from a static pipe flow case and then continuing with an impulsive
injection of the fluid through the nozzle, we first obtain solutions for the
velocity and polymeric stress fields and as the flow approaches steady state
we compare the radial profiles of the dimensionless axial velocity field, and
axial and shear polymeric stresses as they evolve along the nozzle to the
analytical solutions for these fields, which have been derived assuming
steady-state and fully-developed axisymmetric flow (Tomé et al., 2007; Yapici
et al., 2009). As we show in Figure 2 (a, b); we obtain excellent agreement
between the computed velocity and stress fields and the analytic results.

Table 1: Parameter values utilised to validate the flow predictions within the nozzle using the immersed boundary method. $De$ | $Oh$ | $\beta$ | $L^{2}$ | $We$ | $\ell_{nozzle}$ | $LVL$
---|---|---|---|---|---|---
1 | 0.2 | 0.85 | $\infty$ | 16 | 4 | 14

(a) (b)

Figure 2: Validation of the predicted fully-developed flow profiles within the
nozzle with the parameter values given by Table 1. The analytical solution and
corresponding simulation data for the dimensionless axial velocity component
and the polymeric stress components are shown in (a) and (b), respectively, as
a function of dimensionless position across the nozzle. Here, for the
analytical solution (Tomé et al., 2007; Yapici et al., 2009) the dimensionless
strain rate is $Wi=\tau U_{0}/R_{0}=4$.

Following the validation of the steady pipe flow case, we now investigate the
flow inside the nozzle when an axial pressure oscillation with an amplitude
$\epsilon_{p}=0.4$ is forced at the inlet using the expression given by Eq.
(14) with a dimensionless wavenumber $k=0.6$ (with tildes suppressed):

$-{\boldsymbol{\nabla}}{p}\cdot\mathbf{\hat{n}_{I}}=8Ca\left(1+\epsilon_{p}\sin\left(\sqrt{We}\
k\ {t}\right)\right),$ (14)

where $\mathbf{\hat{n}_{I}}$ is the unit normal vector to the inlet boundary,
$Ca=\eta_{l}U_{0}/\gamma$ is the capillary number, and
$We=\rho_{l}U_{0}^{2}R_{0}/\gamma$ is the Weber number, which is related to
the Ohnesorge number defined in Eq. (6) through the relation $Ca=\sqrt{We}\
Oh$. The value of $k=0.6$ is selected according to the linear stability
analysis presented in Section 3.2 in order to correspond to the most unstable
growth rate of the imposed perturbation.

In Figure 3(a) we compare the temporal evolution of the radially-averaged
magnitudes of the dimensionless isotropic pressure
$\bar{p}=\int^{1}_{0}p(r,x=0)rdr$ and axial velocity
$\bar{u}_{x}=\int^{1}_{0}u_{x}(r,x=0)rdr$ at a position close to the nozzle
inlet. Specifically, Figure 3(a) indicates the establishment of a phase lag
and a corresponding amplitude deviation between the pressure and axial
velocity component, as expected from the analysis of (Womersley, 1955). In
Figure 3(b) we also show the evolution of the radially-averaged dimensionless
elastic stress components
$\bar{\sigma}_{p,xx}=\int^{1}_{0}\sigma_{p,xx}(r,x=0)rdr$ and
$\bar{\sigma}_{p,rx}=\int^{1}_{0}\sigma_{p,rx}(r,x=0)rdr$, validating in both
cases the smooth periodic evolution of these quantities. In particular, we
show that initially both the polymer shear stress and streamwise axial stress
evolve together after the flow is impulsively started at short times $t\leq
5$. However, as time increases, the radially-averaged axial and shear
polymeric stress components in the nozzle exhibit distinct trends in magnitude
with the dimensionless axial stress significantly increasing and at long times
($t\geq 5$), dominating the shear stress.

(a)

(b)

Figure 3: Temporal evolution of (a) the cross-sectionally averaged axial velocity component $\bar{u}_{x}$ and pressure $\bar{p}$; and (b) cross-sectionally averaged axial $\bar{\sigma}_{p,xx}$ and shear $\bar{\sigma}_{p,rx}$ polymeric stress components in the nozzle. The parameter values are given in Table 1 and the pressure gradient is set by Eq. (14) with $\epsilon_{p}=0.4$ and $k=0.6$. Table 2: Simulation parameters for the jetting process. The ratios $\rho_{g}/\rho_{l}$ and $\eta_{g}$/$\eta_{l}$ are the same as those used by Turkoz et al. (2018). The dimensionless wavenumber $k$ is selected to be the value that corresponds to the most unstable mode of perturbation expected from linear stability analysis.  $De$ | $Oh$ | $\beta$ | $L^{2}$ | $We$ | $\rho_{g}/\rho_{l}$ |  $\eta_{g}$/$\eta_{l}$ | $k$ | $\epsilon_{p}$ | $\ell_{nozzle}$ | $\ell_{domain}$ | $LVL$
---|---|---|---|---|---|---|---|---|---|---|---
1 | 0.2 | 0.85 | 900 | 16 | 0.01 | 0.01 | 0.6 | 0.4 | 4 | 100 | 14

(a)

(b)

Figure 4: Simulation results at $t=47.8$ for the parameter values listed in
Table 2 showing (a) the volume fraction with the mesh highlighting the
adaptive mesh refinement in the near-interface regions. The interface is shown
with a light green line and the mesh attains its maximum density in the
regions which contain thin fluid threads. The spatial distributions of the
axial velocity and polymeric stress components are shown in the following
contour plots; (b) Local enlargements in the range of $34\leq x<41$ of the
spatial distributions of the volume fraction with the mesh, the axial
velocity, and polymeric stress components in the polymeric thread highlighting
the highly localised distribution of the elastic stress and fore-aft asymmetry
of the satellite drop that develops far downstream of the nozzle.

##  3 Results and Discussion

###  3.1 Jet evolution and breakup

We present numerical simulations of a low-speed, axisymmetric viscoelastic jet
with $De=1$, $Oh=0.2$, $\beta=0.85$, and $We=16$, including the flow within
the nozzle, where the mesh resolution is gradually increased in the range
$12<LVL<14$ $\left(0.02>\Delta x_{minimum}>0.006\right)$; the rest of the
parameters are given in Table 2. Figure 4(a) shows a contour plot of the
volume fraction of fluid in the domain at $t=47.8$, as well as the spatial
distribution within the fluid phase of the dimensionless axial components of
the velocity, $u_{x}(r,x)$, and polymeric stress field, $\sigma_{p,xx}(r,x)$.
The capillarity-driven deformation of the jet is evident as it exits the
nozzle, and this leads ultimately to drop formation. The contour plot of the
axial velocity field demonstrates the initial parabolic profile of the axial
velocity component which develops upstream of the nozzle exit plane, as
expected from Figure 2(a) discussed in section 2. This contour plot also shows
that material elements at the centerline of the liquid jet leave the nozzle at
the maximum velocity and the velocity field rapidly rearranges (within one
perturbation wavelength) so it is then advected downstream at a uniform
average speed. The spatial distribution of the axial polymeric stresses
responds more slowly. Inside the nozzle, the polymeric chains undergo strong
shearing close to the wall of the pipe due to the no-slip boundary condition
resulting in large stresses. Downstream of the nozzle exit plane, a zero shear
stress interfacial condition is imposed, which replaces the no-slip condition
on the inside of the nozzle walls. The axial elastic stress component relaxes
within the beads that form as the perturbed interface evolves under the action
of capillarity, but locally increases in the thin ligaments that develop. This
local increase is driven by the large capillary pressure in the filament as it
thins towards breakup.

Figure 4(b) offers a closer inspection of the simulation results in Figure
4(a). In particular, Figure 4(b) reveals that the two leading drops are
separated by a thin filament on which a much smaller, satellite drop has
formed. These “beads-on-a-string structures” (BOAS) are the result of a
balance between capillarity and the viscoelasticity of the polymer and have no
direct analogue in low-speed jets of Newtonian fluids undergoing deformation
and breakup. It is also clear that the AMR scheme within _Basilisk_ has been
deployed appropriately for refinement close to the interface and to resolve
accurately the stresses in these thin string-like filaments. We also note that
the thread is not perfectly fore-aft symmetric: as time evolves the left upper
side of the thread (i.e. the side which is closest to the nozzle) is observed
to move slightly faster than the lower side (furthest from the nozzle), and
the satellite droplet moves downstream towards the leading drop. Additionally,
the dimensionless axial polymeric stresses are seen to attain very high values
up to $\text{max}(\sigma_{p,xx})\approx 50\gamma/R_{0}$ in the thin ligaments
that develop on either side of the smaller droplet due to the response of the
polymer molecules to the elongational flow. Inside the satellite drop and the
two bigger drops, the stresses relax to values close to zero.

(a)

(b)

Figure 5: (a) A sequence of jet profiles $R(x,t_{i})$ highlighting the
formation of beads-on-a-string structures, as well as primary and satellite
drops, captured each period when the velocity attains its minimum values in
the inlet (these times denoted $t_{i}$ can be seen in Figure 3(a)). (b) Space-
time diagram or “kymograph” showing $R(x,t)$ as each wave pulse is ejected
from the nozzle, leading to the formation of a new primary bead as it flows
away from the nozzle. The simulation parameters are provided in Table 2.

Figure 5(a) presents an alternative Lagrangian view of the jetting process. We
show the temporal evolution of the interface plotted at a sequence of times
denoted $t_{i}$ when the oscillating velocity forcing of the injection at the
inlet attains its minimum value (i.e. $t_{i}\approx 2.3+2\pi n/(\sqrt{We}\ k)\
\text{with}\ n=0,1,2,...$). Mesh resolution represented by different $LVL$
values increases as time increases to ensure that the dynamics are captured
accurately. Each pulsed wave-like disturbance that emanates from the nozzle
results in a local necked region (observable by the blue contours in Figure
5(b)) associated with each new emerging primary bead that travels downstream
at a constant velocity. After an initial transient period of developing flow,
periodicity is observed in terms of the locations where a droplet is formed,
as well as where the thin fluid ligament connects the leading droplet to the
rest of the jet. Figure 5(b) shows the corresponding “kymograph”, which
highlights the periodicity over the entire spatio-temporal spectrum of the
thinning dynamics of the viscoelastic jet. In particular, the contour plot
demonstrates how the local minima in the radius of the jet evolve both in
space and time, where the colour scale ranges from the initial value of the
radius, $R_{0}$, down to the minimum cell size ($\approx 0.6\%R_{0}$ for
$LVL=14$). The kymograph also permits us to track the formation and the
detachment of the leading droplet, observed at early times after which
periodicity of the jet evolution is established. The spatiotemporal
development of satellite droplets, which are represented by cyan-coloured
streaks of smaller radius and travel downstream at almost constant speed is
also highlighted, while the dark blue regions in the contour plot demonstrate
the complete detachment of each of the formed droplets.

Figure 6: Temporal evolution of the jet radius associated with four distinct
regions: ‘I’ and ‘II’ close to the nozzle and at intermediate distances from
it wherein the dynamics are linear and weakly nonlinear, respectively. Regions
‘III’ and ‘IV’ are further away from the nozzle where the dynamics are
strongly nonlinear featuring the formation of the first satellite beads and
beads-on-a-string structures, respectively. The contour plot depicting the
shape of the jet at $t=47.8$ in regions ‘I’-‘IV’ shows the spatial evolution
of the axial component of the polymeric stress tensor accompanied by the
corresponding scale of the $x$-axis. The parameter values are the same as in
Table 2.

In Figure 6, we show the temporal evolution of the jet radius at four fixed
Eulerian locations along the jet axis, which correspond to four distinct
regions of the breakup process (labelled ‘I’ - ‘IV’) and are characterised by
their proximity to the nozzle exit plane and the nature of the evolution of
$R(x,t)$. The corresponding jet profile is coloured by the magnitude of the
axial component of the polymeric stress tensor. Relatively close to the exit
plane of the nozzle (which is located at $x=4$) in region ‘I’ ($4\leq x<14$)
the jet radius exhibits essentially linear dynamics (Middleman, 1965; Goldin
et al., 1969; Brenn et al., 2000) characterised by a sinusoidal response to
the pressure gradient forcing set by Eq. (14).

Following this linear phase, Figure 6 shows that in a second region denoted
region ‘II’ ($14\leq x<24$) the jet radius exhibits a weakly nonlinear
behaviour. Further away from the nozzle ($24\leq x<34$) the radial
perturbation strongly enters into a nonlinear regime denoted region ‘III’
featuring bead formation separated by thin filaments, as shown in Figure 6.
The dynamical evolution of $R(x,t)$ at $x=35$ in region ‘IV’ ($34\leq x<41$)
are strongly nonlinear; here, in the elasto-capillary regime, the thin
ligaments are deformed by capillary forces while elastic stresses delay
interface breakup through the development of stable viscoelastic threads and
beads-on-a-string structures. It is also clear from a close inspection of
Figure 6 that the initial transient response (corresponding to when the first
“leading” droplet exits the nozzle and then passes each Eulerian location)
takes longer to decay further from the nozzle but eventually the jet attains a
perfectly periodic structure (corresponding to the diagonal lines observed in
the space-time diagram of Figure 5).

In what follows, we first use linear stability theory to study the regions
closest to the nozzle, i.e. region ‘I’ and ‘II’, before embarking on a
detailed analysis of regions, ‘III’ and ‘IV’, in which the dynamics become
increasingly nonlinear. In experimental visualisation of the jet breakup, it
is common to follow the evolution of local maxima in the jet radius (leading
to the formation of primary beads) as well as local minima in the necks (which
ultimately lead to pinch-off). The situation is more complex in viscoelastic
jets because the large elastic stresses that develop in the neck can inhibit
or totally prevent pinch-off. It is thus important to always follow the same
Lagrangian element when attempting to relate local rates of thinning to
material properties such as the local extensional viscosity in a material
element. To aid our analysis of the nonlinear dynamics, we move to a
Lagrangian description so that we follow a specific material element as it is
ejected and transported away from the nozzle and is increasingly deformed by
capillary effects in space and time. We show in Figure 7 the space-time
evolution of the jet interface (oriented horizontally here to conserve space)
highlighting the evolution of the wave crests and troughs with distance from
the nozzle. This representation allows to ensure that the same minimum is
followed in space and time, which is essential for estimating the local
instantaneous rate of thinning in region ‘IV’ accurately. As the jet thins and
the beads-on-a-string structure fully develops (at $x\simeq 35$), two thin
threads are formed, one on each side of the satellite droplet. This results in
two local minima in the jet radius $R_{min1}^{[\alpha]}$ and
$R_{min2}^{[\alpha]}$ following each labelled primary bead $[\alpha]=$
A,B,C,…, for each period of the upstream forcing. These two emerging minima
are henceforth labelled ‘min1’ and ‘min2’ for each local minimum in every neck
region.

Figure 7: Space-time plot of the spatio-temporal evolution of the jet
highlighting the detection and tracking of two different minima of the jet
radius on either side of a secondary droplet, always observed for $x\geq 25$.
Each local maximum corresponding to the formation of a primary bead is shown
by a filled black circle symbol with a distinctly coloured border, as well as
a Lagrangian label $[\alpha]$. Each local minimum in the thinning fluid
ligament between neighbouring primary beads is shown by a filled
(‘$R_{min1}^{[\alpha]}$ \- min1’) and hollow (‘$R_{min2}^{[\alpha]}$ \- min2’)
symbol of the corresponding colour for each neck region established behind
each primary bead $[\alpha]$. The enlarged view shows the jet profile between
two beads with Lagrangian labels $[\alpha]=$ N, O. The velocity at which the
viscoelastic jet evolves downstream is also indicated. The simulation
parameter values are the same as in Table 2.

(a)

(b)

Figure 8: Comparison of the numerical predictions with those from a linear
stability analysis for regions ‘I’ and ‘II’: (a) semi-log plot of the temporal
evolution of the deviation of the jet radius from its base state for an
imposed perturbation given by Eq. (15) with $\delta=0.13$ and $k=0.6$. The
rest of the simulation parameters remain unchanged from Table 2. Here,
$t-t_{p}$ corresponds to the elapsed time since a specific local maximum in
the jet radius which coincides with the formation of a primary bead labelled
$[\alpha]=$ L, M, N, O, P is formed as a fluid ligament flows away from the
nozzle. Regions ‘III’ and ‘IV’ are characterised by nonlinear dynamics for
which linear theory is not appropriate; regions ‘I’-‘IV’ map onto those
identified in Figure 6. (b) Dispersion curves generated via the solution of
Eq. (17) (Middleman, 1965; Brenn et al., 2000) for various Weber numbers and
wavenumbers; simulation data for the specific local neck regions highlighted
in Figure 7 were used for the filled circle symbol at $k=0.6$ and the
simulation parameters provided in Table 2; the hollow symbols have been
generated for the same $Oh,De$ and $We$ ($Oh=0.2,De=1,We=16$) conditions but
different $k$ values.

###  3.2 Region ‘I’ and ‘II’: Linear and weakly nonlinear evolution of jet
profile

In regions ‘I’ and ‘II’, we use linear stability theory to characterise the
thinning dynamics. Here, the decrease in the dimensionless jet radius is
expected to be described by a growing perturbation of the general form:

$R_{min}(t)=1-\delta\exp(\Omega t),for\ t\geq 0$ (15)

where $\delta$ is the initial perturbation amplitude and
$\Omega=\Omega_{r}+i\Omega_{i}$ is the complex growth rate which depends
parametrically on $Oh$, $De$, and $k$: $\Omega_{r}>0$ indicates the presence
of linear instability. In Figure 8(a) we show a semi-log plot of
$1-R_{min}^{[\alpha]}$ as a function of time whence we have subtracted an
interval $t_{p}$ that represents the instant when each of the identified wave
pulses, which lead to the formation of primary beads labelled $[\alpha]=$ L,
M, N, O, P, exited the nozzle. We also illustrate the locations of the
stationary Eulerian points located at 5, 15, 25, and 35 nozzle diameters from
the injection inlet. Given the velocity of the jet, each of these fixed points
can be associated with a specific value of $t-t_{p}$ which corresponds to the
time when a specific material element passes through one of each of these
locations. Inspection of this plot reveals that even though the perturbations
in region ‘I’ are small as illustrated in Figure 6, the residual stresses in
the jet, the rearrangement of the velocity profile in the jet (from parabolic
to plug-like) and the pinning conditions of the free surface to the nozzle
exit at $x=4$ all influence the local growth rate of perturbations. Hence, it
is the second region denoted ‘II’ which is best characterised by the linear
theory.

In region ‘II’ (beyond approximately one jet diameter from the nozzle), the
exit boundary conditions have been forgotten and small perturbations to the
radius grow exponentially under the action of capillary squeezing. Figure
(8)(a) shows that region ‘II’ remains linear (on a semi-log plot) over a
sufficiently large time interval that it is possible to calculate the gradient
which can be compared with the dimensionless growth rate $\Omega_{r}$ in Eq.
(15). The growth rates obtained from the numerical simulation for
$Oh=0.2,De=1$ and $We=16$ and linear stability theory (Brenn et al., 2000)
(see Appendix B for details) are $\Omega_{r}\approx 0.246\pm 0.006$ and 0.25,
respectively, for $k=0.6$, demonstrating excellent agreement, as is also shown
in Figure 8(b).

In particular, we can compare the numerically predicted growth rate with those
obtained from linear theory for a range of $k$ values. Figure 8(b) shows the
resulting dispersion curves computed using the analysis of Middleman (1965)
for a viscoelastic liquid jet in an inviscid gaseous environment.
Specifically, the dispersion curve of an inviscid Newtonian jet ($Oh=0,De=0$)
with no inertia ($We=0$) (Rayleigh, 1879) is first presented (blue dashed-
dotted line), while viscosity effects are then added for the Newtonian jet
($Oh=0.2,De=0,We=0$) resulting in a strongly reduced growth rate of the
perturbation (orange dashed line). Subsequently, the effects of
viscoelasticity ($De=1$) are considered (solid lines): first, the influence of
inertia is neglected ($Oh=0.2,De=1,We=0$) showing that the linear viscoelastic
jet is more unstable (slightly larger maximum growth rate) than a Newtonian
fluid of the same viscosity. That is, elasticity contributes to a little
faster rate of thinning in the neck of the liquid jet in the linear regime
(Middleman, 1965). Brenn et al. (2000) expanded the work of Middleman (1965)
by incorporating the effects of the momentum flux arising from the injection
of the jet (ie Weber numbers $We>0$. All of the corresponding dispersion
curves ($Oh=0.2,De=1,We>0$) exhibit positive growth rates at $k=1$ instead of
$\Omega_{r}=0$ as obtained when $We=0$. In addition, the presence of fluid
inertia is destabilising, leading to a higher maximum wave growth rate, a
shift of the most unstable mode to larger $k$ values, and a wider range of
wavenumbers for which $\Omega_{r}>0$.

Comparison of the linear theory predictions to the computed growth rates
resulting from the numerical simulations for $Oh=0.2$, $De=1$, $We=16$ and
$k=0.3,0.4,0.5,0.6,0.7$ and $0.8$ (obtained using a similar procedure to that
discussed in connection to Figure 8(a)), shows good agreement. Specifically,
all the circular symbols in Figure 8(b) result from the analysis of the
corresponding rates of change, equivalent to the analysis shown in Figure
8(a). Analysis of at least five different Lagrangian wave pulses exiting the
nozzle result in the error bars shown in the figure. Therefore, the slight
deviations seen in Figure 8(b) between the simulation predictions and the
dispersion curve computed for $Oh=0.2,De=1$ and $We=16$ (green curve) indicate
the limits of our ability to resolve small differences in the resulting growth
rates due to the small neck-to-neck variations in the profiles of the thinning
jet considered in each of the simulations.

(a)

(b)

Figure 9: (a) Thinning dynamics and (b) the temporal evolution of the local
Weissenberg number, $Wi$ in each Lagrangian thinning neck region. Data from
each of the local neck regions depicted in Figure 7 are used (for simulation
parameters listed in Table 2). In (b), the plateau values corresponding to
$Wi=2/3$ ($\dot{\epsilon}_{min}=2/3\tau$) and $Wi=1$
($\dot{\epsilon}_{min}=1/\tau$), which coincide with the two distinct limits
to the thinning dynamics in the elasto-capillary regime (Keshavarz et al.,
2015; Mathues et al., 2018), are also shown. The regions ‘I’-‘IV’ identified
in (a) map directly onto those discussed in Figures 6 and 8(a). Points
labelled ‘P1’-‘P4’ in (b) correspond to the ends of regions ‘I’-‘IV’ in (a),
respectively.

###  3.3 Regions ‘III’ and ‘IV’: Nonlinear evolution of jet profile

The evolution in jet profile in the nonlinear regimes is first analysed by
tracking the temporal decrease of the jet radius following local minima in
each neck region between primary beads in a Lagrangian way as it travels
downstream away from the nozzle, as shown in Figure 9(a), for the same
simulation parameters as in Table 2. This highlights the emergence of four
distinct regimes, as defined in Figure 6: ‘I’, ‘II’, ‘III’, and ‘IV’
characterised by $R_{min}\geq 0.75$, $0.75>R_{min}\geq 0.25$,
$0.25>R_{min}\geq 0.1$ and $0.1>R_{min}\geq 0.006$, respectively, where
$0.6\%R_{0}$ is the mesh resolution limit according to the maximum $LVL$ value
achieved in this work. As the thinning jet enters region ‘IV’ it is clear that
the radius of the local minima ‘min1’ and ‘min2’ (corresponding to the local
minimum jet radius in each of the two thin ligaments between beads) decreases
exponentially in time. This is the elasto-capillary (EC) regime anticipated in
the Introduction. Finally, at very small radii below $R_{min}\leq 0.04$, there
is a deviation from the exponential thinning corresponding to the onset of
finite extensibility effects. In this regime, the thread radius is ultimately
expected to thin linearly in time resulting in finite time breakup (Entov &
Hinch, 1997; Renardy & Renardy, 2004).

To determine the characteristic time-scale that best describes the exponential
thinning, the temporal evolution of the local dimensionless strain rate $Wi$,
defined in Section 2.1, is plotted in Figure 9(b). This plot is characterised
by four points, ‘P1’-‘P4’, highlighting the non-monotonic evolution of the
local strain rate in the thinning jet in agreement with what has been shown by
Tirtaatmadja et al. (2006) and more recently by Rajesh et al. (2022). The
points ‘P1’-‘P4’ are respectively associated with the ends of regions
‘I’-‘IV’, identified in Figure 9(a): ‘P1’ is associated with the end of region
‘I’ in which the influence of the exit nozzle on the thinning jet radius is
felt, while ‘P2’ coincides with the end of region ‘II’ and exponential growth
in disturbances evident in Figure 8(a). In the transition regime (region
‘III’), the local dimensionless strain rate $Wi$ associated with the evolution
in the local minima ‘min1’ and ‘min2’ passes through a local maximum and then
decreases until Point ‘P3’, which corresponds to the end of region ‘III’.
Point ‘P3’ is characterised by a local value of the Weissenberg number, which
we shall generically denote by $Wi_{EC}$, that remains approximately constant
for a period of time whose duration depends on We and $L^{2}$, as will be
discussed below, and heralds the transition to the elasto-capillary regime in
region ‘IV’. Towards the end of region ‘IV’ the polymer molecules reach their
maximum extensibility and the local Weissenberg number undergoes a steady
increase. In this regime, beads-on-a-string structures are also formed
reflecting the complex balance between capillary and elastic stresses. For the
relatively large We and small $L^{2}$ values used to generate Figure 9(b) we
see that the plateau value of $Wi_{EC}\approx 1$ is only established for a
narrow time interval and then the strain rate steadily increases until point
’P4’, which coincides with the end of region ‘IV’ when the finite
extensibility limit is reached and the local strain rate diverges as the local
radius of the thread decreases to zero.

We note that data from the minimum radius associated with the necks
established behind five different Lagrangian primary beads (labelled L
$\rightarrow$ P) are shown in Figure 9(a,b). It is clear that the jet exhibits
self-similar thinning dynamics in the initial inertio-capillary regime; this
is evidenced by the overlapping curves in Figure 9(a) in regions ‘I’ and ‘II’.
If nonlinear elastic effects were not important then the linear perturbations
to the jet radius would continue to grow exponentially until the radius
locally approaches zero (according to Eq. (15)) and the local strain rate
$\dot{\epsilon}_{min}=-(2/R_{min})(dR_{min}/dt)$ would evolve as indicated by
the curves labelled ‘Linear Stability’ in Figures 9(a,b). However, as
nonlinear viscoelastic effects in the fluid thread become important, the rate
of thinning decreases and the local strain rate passes through a maximum
labelled shortly after point ’P2’. The local thread dynamics are also self-
similar in both the transition region and the elasto-capillary regime as
demonstrated by the superposition of curves in regions ‘III’ and ‘IV’, as well
as ‘P2’ and ‘P3’ points in Figures 9(a,b), respectively. It is also
interesting to note that the local maximum in the strain rate associated with
the minimum radius denoted ‘min1’ (the leading ligament ahead of the forming
satellite bead) is higher than the one associated with ‘min2’, thereby
highlighting the role of the momentum flux in the jet.

Figures 10(a,b) show the dynamical evolution of the axial velocity along the
centerline and the polymeric stress component at the jet centerline for two
Lagrangian points corresponding to ‘min1’ and ‘min2’. Also shown as inserts
are snapshots of the axial velocity contours and polymeric stress fields taken
at times that correspond to regions ‘I’, ‘II’, ‘III’, and ‘IV’, as defined in
Figures 6 and 9(a).

(a)

(b)

Figure 10: Temporal variation of (a) the axial velocity component, and (b) the
axial component of the elastic stress. Both are evaluated at the jet
centreline using data from the local minima associated with the necks as
depicted in Figure 7 considering the same simulation parameters provided in
Table 2. Also shown in (a) and (b) are profiles which depict the shape of the
jet, coloured by contours indicating the magnitude of the axial components of
the velocity and the elastic stress, respectively, during times:
$0<t-t_{p}\leq 1$, $1<t-t_{p}\leq 4.4$, $4.4<t-t_{p}\leq 6.2$ and
$6.2<t-t_{p}\leq 7.8$ that correspond to regions ‘I’-‘IV’ presented in Figure
9(a), respectively.

The centerline velocity associated with both ‘min1’ and ‘min2’, which
correspond to the same Lagrangian element in the initial stages (regions
‘I’-‘II’), first undergoes a decrease in region ‘I’ as the jet exits the
nozzle and the velocity profile rearranges, followed by a slow increase in
region ‘II’ associated with local perturbations growing according to the
linear stability analysis. In the transition regime, region ‘III’, the two
local minima in the thread radius behind each primary bead start following
different dynamics as the beads-on-a-string structure starts forming: the
axial velocity associated with point ‘min1’ (‘min2’) decreases (increases)
until reaching a minimum (maximum). Shortly after reaching their local
extremal value, the axial velocity in Figure 9(a) then exhibits an
approximately linear decrease (increase) during region ‘IV’ as each part of
the thinning jet approaches a constant velocity.

In contrast, the axial component of the polymeric stress shown on a semi-log
scale in Figure 10(b) exhibits a more complex and non-uniform rise over time
with high rates of increase at the exit plane of the nozzle in region ‘I’, and
in region ‘III’ while a slower increase occurs in the region of linear
instability (region ‘II’). In the elasto-capillary regime (region ‘IV’), the
exponential increase in tensile stress within the thinning filament over time
is clearly seen. As expected this results in the largest values of the axial
component of the polymeric stress in the thin ligament (Clasen et al., 2006;
Eggers et al., 2020; Deblais et al., 2020). Finally beyond ‘P4’ the finite
extensibility limit is approached and the radius decreases to zero as the
local strain rate and resulting stress in the thinning thread diverge. Once
again, we note that all the data obtained from the five individual Lagrangian
local necked regions used to generate Figure 10 collapse to form a master
curve for the evolving axial velocity and stress components at the jet
centreline highlighting the periodicity and self-similarity of the established
dynamics.

In contrast to the free viscoelastic filament undergoing thinning in the
absence of a mean flow (the $We=0$ case) presented recently in Turkoz et al.
(2018), in the jetting process for a low to moderate Weber number, the
momentum flux of the ejected fluid at the nozzle exit stimulates the
exponential decrease of the jet radius and the development of a fore-aft
asymmetric beads-on-a-string structure, as indicated by the distinct evolution
of the two local minima that separate the formation of the viscoelastic
ligament in Figure 10(a). Below we further investigate the role of the
injection flowrate in the pinch-off dynamics and the resulting relative
contributions of inertia, capillarity, and elastic stresses in the evolving
viscoelastic thread.

###  3.4 Effect of Inertia

In Figure 11 we show the effect of altering the injection rate on the
resulting jet dynamics by varying the Weber number; we plot snapshots of the
interface shape coloured by the magnitude of the axial polymeric stress
component field for $We=8$, $16$, and $36$ with the rest of the parameters
remaining unchanged from those shown in Table 2. The smallest Weber number
studied was chosen to be larger than that associated with the so-called
“gobbling limit” (typically seen around $We\approx 2$ (Clasen et al., 2009)).

Figure 11: The effect of fluid inertia on the interfacial thinning dynamics:
contour plots of the jet shape coloured by the magnitude of the axial
component of the elastic stress for $We=8$, 16, and 36 at $t=53.6$, 47.8, and
29.2 respectively, with the rest of the parameters remaining unchanged from
Table 2. Also shown are enlarged views of the leading bead and the
interconnecting ligament regions for each value of Weber number.

It is seen that the jet length increases with We and the thinning dynamics are
accompanied by a concomitant rise in the number of undulations that develop
into necks with longer strings separating the formed beads. According to the
dispersion curves in Figure 8(b), the magnitude of the Weber number has a weak
influence on the instability growth rates, and the perturbations are therefore
advected further away from the nozzle before entering the elasto-capillary
regime when the Weber number is increased. Furthermore, the size of the
satellite drops along the ligaments interconnecting the primary drops also
increases with We and their position is shifted downstream towards the leading
bead; the size of the primary beads, however, only appears to be weakly
dependent on We. Moreover, from the contour plots of the elastic stresses
shown in Figure 11, it is clear that the increase in We results in higher
polymeric stresses within the nozzle and correspondingly at the exit, but
these largely relax within a few jet diameters and there is only a slight
increase in the stress levels attained in the thin viscoelastic threads.

We also study the temporal evolution in the local dimensionless strain rate
$Wi(t)=\tau\dot{\epsilon}_{min}(t)$ in Figure 12(a) for Weber numbers 8, 16,
and 36 with the rest of the parameters remaining unchanged from Figure 10. As
the profiles in Figure 10 are identical for each neck established behind a
formed primary bead labelled L, M, N…, we focus only on one local neck
henceforth in Figure 12. We also consider the flow dynamics only after
dimensionless times ($t-t_{p}\geq 1$) during which the effect of the nozzle
exit becomes less pronounced. In each case, it is clear that the evolution in
the local strain rate in a fluid neck, as it evolves along the jet, shows all
the features documented in Figure 9 with a slow increase in Wi(t) as the
disturbances grow, a local maximum in the deformation rate before the necking
material element enters the elasto-capillary (EC) regime (region ‘IV’) in
which the Weissenberg number approaches a locally constant value that we
denote $Wi_{EC}$. However, it can be seen that decreasing the level of inertia
in the jet results in the approach to a plateau value $Wi_{EC}=2/3$, in marked
contrast to the $We=16$ and $36$ cases where $Wi_{EC}=1$. It is also observed
that increasing the Weber number leads to larger values of the local maxima in
the Weissenberg number obtained after point ‘P2’, which coincides with the
transition to the elasto-capillary regime. It is also clear that in the case
of the smaller Weber number ($We=8$) that the transition from the
characteristic points labelled ‘P3’ to ‘P4’ in Figure 10 is more gradual
compared to $We=16$ and $36$.

(a)

(b)

Figure 12: The effect of increasing fluid inertia on (a) the temporal
evolution of the local dimensionless strain-rate
$Wi(t)=\tau\dot{\epsilon}_{min}(t)$ in the local thinning neck of the fluid
jet, and (b) the capillary, inertial, and elastic contributions to the
simplified one-dimensional force balance (Eq. (16)), for $We=8$, 16, and 36
with the rest of the parameters remaining unchanged from Table 2.

In order to assess the effect of jet inertia on the relative contributions of
the capillary, viscous, inertial, and elastic stresses, we consider a simple
one-dimensional model for dimensionless axial force along the jet, $F_{x}$,
given by Mathues et al. (2018) in which the viscous forces are neglected, and
results in the expression:

$\frac{F_{x}}{\pi
R_{min}^{2}}\sim\frac{1}{R_{min}}+\frac{\dot{R}_{min}^{2}}{2}+\bar{\sigma}_{p,xx},$
(16)

where the terms on the right-hand side correspond to the dimensionless
capillary pressure, inertia, and average axial polymeric stress in the radial
direction.

In Figure 12(b), we plot the contributions to $F_{x}/\pi R_{min}^{2}$ obtained
from each of the terms on the right-hand side of Eq. (16). In each case
initially, capillary pressure dominates while the initial magnitude of the
average elastic stresses is high due to the significant effects of the nozzle.
Then, as the boundary condition at the exit of the nozzle changes from no-slip
to free-slip, the average axial polymeric stress relaxes in all cases, before
it grows as the neck thins and ultimately (after a period of exponential
growth) dominates the axial force along the jet. Inertial terms are always
small but become increasingly unimportant as the jet and length scales become
smaller. Increasing We enhances slightly the inertial contribution as we
approach thread pinch-off but with the balance of the capillary action in the
thread by the elastic effects being established significantly more quickly.
The level of elastic tensile stresses that can be sustained in the thinning
fluid thread also depends on the finite extensibility of the dissolved polymer
molecules. In the next section, we examine the effect of changing the value of
the FENE parameter $L^{2}$ on the temporal evolution of the local stretching
rate and the magnitude of the tensile stress difference that develops in the
thinning fluid ligaments to assess its influence on the jet behaviour in the
elasto-capillary regime and the resulting jet lifespan.

###  3.5 Effect of $L^{2}$

(a)

(b)

Figure 13: The effect of increasing extensibility parameter $L^{2}$ on (a) the
temporal evolution of the local dimensionless strain rate
$Wi(t)=\tau\dot{\epsilon}_{min}(t)$, and (b) evolution of the capillary,
inertial, and elastic contributions to the simplified force balance (Eq. (16))
for $L^{2}=900$, 2500, 10000, and $We=16$, with the rest of the parameters
remaining unchanged from Table 2. Also shown for comparison in (a) the
Oldroyd-B limit corresponding to $L^{2}\rightarrow\infty$.

In Figure 13(a) we show how the dimensionless strain rate fin a representative
fluid neck varies with the extensibility parameter $L^{2}$ at a Weber number
$We=16$ (the rest of the parameters remain unchanged from Table 2). It is
clear that an increase in the extensibility of the polymeric chains beyond
$L^{2}=900$ leads to the strain rate in the elasto-capillary regime converging
progressively to a plateau of value $Wi=2/3$ for a time duration that appears
to be weakly dependent on $L^{2}$. In contrast, for $L^{2}=900$, as discussed
above (see Figure 9(b)), the limited extensibility of the chains prevents a
full elasto-capillary balance from being established and there is a rapid
divergence in the local strain rate from point ‘P3’ towards ‘P4’, with the
$Wi=2/3$ plateau never being approached. Moreover, the local peaks in $Wi$,
which coincide with the transition to the elasto-capillary balance, increase
with $L^{2}$, reaching saturation as the Hookean dumbbell limit
$L^{2}\rightarrow\infty$ is approached. In Figure 13(b) we show that for
$L^{2}>900$, a broader region ‘IV’ (characterised by an exponential increase
of the axial polymeric stress component arises) with a lower rate of change
than in the $L^{2}=900$ case. This shows a more stable local balance between
capillarity and elastic stresses as the axial polymeric stress component can
keep increasing exponentially as a result of the increased extensibility.

In Figure 14 we construct a flow map in $(L^{2},We)$ space in which we collect
the results presented in Figures 12 and 13. The map is coloured by the
magnitude of $Wi_{EC}$ established during the elasto-capillary balance. The
values of $Wi_{EC}$ are computed from the strain rate in the necking filament
at the onset of the elasto-capillary regime and serve to highlight whether or
not the thinning dynamics are significantly accelerated beyond the value
$Wi_{EC}=2/3$ expected in the classic elasto-capillary balance (Entov & Hinch,
1997) depending on We and $L^{2}$.

Figure 14: Flow regime map in $(L^{2},We)$ parameter space coloured by the
magnitude of the Weissenberg number at the start of the elasto-capillary
regime, $Wi_{EC}$, which corresponds to point ‘P3’ labelled in Figure 9(b).
The four contour plots of the jet profile which highlight the shape of the jet
for low and high values of $L^{2}$ and We are coloured by the magnitude of the
axial component of the elastic stress. The rest of the parameters are given in
Table 2.

Additional simulations are performed over a range of We and $L^{2}$ to cover
an extended region of parameter space from low to moderate jet speeds and from
moderate to large polymer chain extensibilities. As indicated by the arrows in
Figure 14, pronounced beads-on-a-string structures are promoted for large We,
associated with longer jet lengths with multiple beads, whilst high values of
$L^{2}$ enable large elastic stresses to develop in the jet and lead to the
formation of longer and thinner ligaments without satellite droplets attached,
as well as slower thinning dynamics. When the axial momentum in the jet is
small and the breakup length of the jet is constrained by the finite
extensibility, an elasto-capillary balance with $Wi_{EC}=2/3$ is established,
similar to the dynamics realised in the CaBER device (Entov & Hinch, 1997;
Anna et al., 2001). However, when the axial momentum in the jet is high and
the length to breakup is very high, the local state of stress in the thinning
neck becomes insensitive to the pinning condition at the exit nozzle. In this
limit, the earlier asymmetric force balance discussed by Clasen et al. (2009)
and Mathues et al. (2018) applies, resulting in a faster local stretching rate
such that $Wi_{EC}\approx 1$.

##  4 Conclusions

We have studied the thinning and breakup of an axisymmetric viscoelastic jet
issuing from a nozzle using the FENE-P constitutive relation which accounts
for finite polymer chain extensibility. We have used the open-source code
Basilisk, which is based on a Volume-of-Fluid interface-capturing methodology
and utilises adaptive mesh refinement for accurate and efficient free surface
flow solutions. The free surface evolution of the jet is coupled to the
upstream flow and the initial polymeric stress development inside the nozzle
by employing a simplified Immersed Boundary Method. The numerical solutions of
the local flow within the nozzle are in excellent agreement with analytical
solutions for the cases of steady and pulsating flows (Womersley, 1955).

Our numerical simulations of the interfacial dynamics capture the development
of beads-on-a-string structures (Clasen et al., 2006) for fixed Ohnesorge and
Deborah numbers over a range of Weber numbers as well as for a range of finite
polymer chain extensibilities, representative of real polymer solutions. We
have highlighted the development of four local regions along the jet axis. The
initial growth of small amplitude perturbations is consistent with linear
stability theory sufficiently close to the nozzle inlet and gives way to
nonlinear dynamics further downstream where the jet evolution is first
governed by the interplay of capillary, inertial, and viscous forces, and
subsequently dominated by an elasto-capillary balance characterised by the
formation of a distinct beads-on-a-string structure. We have also successfully
captured the exponential thinning of the thin highly-stretched ligaments that
form between the primary beads. This elasto-capillary thinning regime is
short-lived for small polymer chain extensibilities but becomes more
pronounced for high $L^{2}$, and the corresponding polymeric tensile stresses
grow larger and larger, in the limit of infinite chain extensibility. In this
limit, adaptive mesh resolution becomes essential and we are able to simulate
the time-dependent evolution of the jet down to minimum feature sizes of
$R_{min}\approx 0.006R_{0}$ and jet lengths as large as $\ell_{jet,max}\approx
100R_{0}$ (see Appendix B for additional details).

Finally, we have explored in detail the local thinning dynamics of the slender
ligaments that develop between the beads-on-a-string structures that evolve
along the jet to resolve differences in previous reports that affect the
determination of a characteristic fluid relaxation time. We construct a flow
map in chain extensibility and Weber number space and calculate the variations
in the local dimensionless extension rate $Wi_{EC}(We,L^{2})$. This map helps
us identify regions of parameter space characterised by the presence or
absence of satellite drop formation, the development of very long viscoelastic
jets with pronounced beads separated by thin strings, and how the thinning
dynamics in the ligaments may vary from $Wi_{EC}=2/3$ (deep red colours) to
$Wi_{EC}=1$ (yellow colour contours). Understanding this systematic evolution
is essential if an accurate value of the characteristic relaxation time in an
unknown fluid is to be extracted from measurements of ligament thinning in a
jetting rheometer or inkjetting device (Morrison & Harlen, 2010; Keshavarz et
al., 2015; Mathues et al., 2018; Xu et al., 2021).

Acknowledgements

This work is supported by the Engineering and Physical Sciences Research
Council, United Kingdom, through the EPSRC PREMIERE (EP/T000414/1) Programme
Grant. Support from Johnson Matthey for Konstantinos Zinelis is also
gratefully acknowledged.

There are no conflicts of interest to report.

## References

  * Amarouchene et al. (2001) Amarouchene, Y., Bonn, D., Meunier, J. & Kellay, H. 2001 Inhibition of the finite-time singularity during droplet fission of a polymeric fluid. Phys. Rev. Lett. 86 (16), 3558–3561.
  * Aniszewski et al. (2020) Aniszewski, W., Saade, Y., Zaleski, S. & Popinet, S. 2020 Planar jet stripping of liquid coatings: Numerical studies. International Journal of Multiphase Flow 132, 103399.
  * Anna et al. (2001) Anna, S. L., McKinley, G. H., Nguyen, D. A., Sridhar, T., Muller, S. J., Huang, J. & James, D. F. 2001 An interlaboratory comparison of measurements from filament-stretching rheometers using common test fluids. J. Rheol. 45 (1), 83–114.
  * Ardekani et al. (2010) Ardekani, A. M., Sharma, V. & McKinley, G. H. 2010 Dynamics of bead formation, filament thinning and breakup in weakly viscoelastic jets. J. Fluid Mech. 665, 46–56.
  * Basaran et al. (2013) Basaran, O. A., Gao, H. & Bhat, P. P. 2013 Nonstandard inkjets. Annu. Rev. Fluid Mech. 45, 85–113.
  * Bazilevesky et al. (1990) Bazilevesky, A.V., Entov, V. N. & Rozhkov, A. N. 1990 Liquid filament microrheometer and some of its applications. Proc. Third Eur. Rheol. Conf. (ed. D. R. Oliver) p. 41 Elsevier Applied Science.
  * Bhat et al. (2010) Bhat, P. P., Appathurai, S., Harris, M. T., Pasquali, M., McKinley, G. H. & Basaran, O. A. 2010 Formation of beads-on-a-string structures during break-up of viscoelastic filaments. Nat. Phys. 6 (8), 625–631.
  * Bhat et al. (2008) Bhat, P. P., Basaran, O. A. & Pasquali, M. 2008 Dynamics of viscoelastic liquid filaments: Low capillary number flows. J. Non-Newt. Fluid Mech. 150 (2-3), 211–225.
  * Bird et al. (1987) Bird, R. B., Curtiss, C. F., Armstrong, R. C. & Hassager, O. 1987 Dynamics of Polymer Liquids Vol. 2. Wiley.
  * Bousfield et al. (1986) Bousfield, D. W., Keunings, R., Marrucci, G. & Denn, M. M. 1986 Nonlinear analysis of the surface tension driven breakup of viscoelastic filaments. J. Non-Newt. Fluid Mech. 21 (1), 79–97.
  * Brenn et al. (2000) Brenn, G., Liu, Z. & Durst, F. 2000 Linear analysis of the temporal instability of axisymmetrical non-Newtonian liquid jets. Int. J. Multiph. Flow 26 (10), 1621–1644.
  * Clasen et al. (2009) Clasen, C., Bico, J., Entov, V. M. & McKinley, G. H. 2009 ’Gobbling drops’: The jetting dripping transition in flows of polymer solutions. J. Fluid Mech. 636, 5–40.
  * Clasen et al. (2006) Clasen, C., Eggers, J., Fontelos, M. A., Li, J. & McKinley, G. H. 2006 The beads-on-string structure of viscoelastic threads. J. Fluid Mech. 556, 283–308.
  * Deblais et al. (2020) Deblais, A., Herrada, M. A., Eggers, J. & Bonn, D. 2020 Self-similarity in the breakup of very dilute viscoelastic solutions. J. Fluid Mech. 904, R2.
  * Dinic & Sharma (2019) Dinic, J. & Sharma, V. 2019 Computational analysis of self-similar capillary-driven thinning and pinch-off dynamics during dripping using the volume-of-fluid method. Phys. Fluids 31 (2), 021211.
  * Dinic et al. (2015) Dinic, J., Zhang, Y., Jimenez, L. N. & Sharma, V. 2015 Extensional Relaxation Times of Dilute, Aqueous Polymer Solutions. ACS Macro Lett. 4 (7), 804–808.
  * Eggers et al. (2020) Eggers, J., Herrada, M. A. & Snoeijer, J. H. 2020 Self-similar breakup of polymeric threads as described by the oldroyd-b model. Journal of Fluid Mechanics 887, A19.
  * Eggers & Villermaux (2008) Eggers, J. & Villermaux, E. 2008 Physics of liquid jets. Reports Prog. Phys. 71 (3).
  * Entov & Hinch (1997) Entov, V. M. & Hinch, E. J. 1997 Effect of a spectrum of relaxation times on the capillary thinning of a filament of elastic liquid. J. Non-Newt. Fluid Mech. 72 (1), 31–53.
  * Entov & Yarin (1984) Entov, V. M. & Yarin, A. L. 1984 Influence of elastic stresses on the capillary breakup of jets of dilute polymer solutions. Fluid Dynamics 19 (1), 21–29.
  * Étienne et al. (2006) Étienne, J., Hinch, E. J. & Li, J. 2006 A Lagrangian-Eulerian approach for the numerical simulation of free-surface flow of a viscoelastic material. J. Non-Newt. Fluid Mech. 136 (2-3), 157–166.
  * Fattal & Kupferman (2005) Fattal, R. & Kupferman, R. 2005 Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation. J. Non-Newt. Fluid Mech. 126 (1), 23–37.
  * Fontelos & Li (2004) Fontelos, M. A. & Li, J. 2004 On the evolution and rupture of filaments in Giesekus and FENE models. J. Non-Newt. Fluid Mech. 118 (1), 1–16.
  * Goldin et al. (1969) Goldin, M., Yerushalmi, J., Pfeffer, R. & Shinnar, R. 1969 Breakup of a laminar capillary jet of a viscoelastic fluid. Journal of Fluid Mechanics 38 (4), 689–711.
  * Keshavarz et al. (2015) Keshavarz, B., Sharma, V., Houze, E. C., Koerner, M. R., Moore, J. R., Cotts, P. M., Threlfall-Holmes, P. & McKinley, G. H. 2015 Studying the effects of elongational properties on atomization of weakly viscoelastic solutions using Rayleigh Ohnesorge Jetting Extensional Rheometry (ROJER). J. Non-Newt. Fluid Mech. 222, 171–189.
  * Lauser et al. (2021) Lauser, K. T., Rueter, A. L. & Calabrese, M. A. 2021 Small-volume extensional rheology of concentrated protein and protein-excipient solutions. Soft Matter 17 (42), 9624–9635.
  * Li & Fontelos (2003) Li, J. & Fontelos, M. A. 2003 Drop dynamics on the beads-on-string structure for viscoelastic jets: A numerical study. Phys. Fluids 15 (4), 922–937.
  * Liu et al. (2023) Liu, L., Guan, X. & Fu, Q. 2023 Numerical simulation study of extensional characteristics impacts on the viscoelastic thread deformation and satellite droplet generation. J. Non-Newt. Fluid Mech. 311, 104955.
  * Lohse (2022) Lohse, D. 2022 Fundamental Fluid Dynamics Challenges in Inkjet Printing. Annu. Rev. Fluid Mech. 54 (1), 349–382.
  * López-Herrera et al. (2019) López-Herrera, J. M., Popinet, S. & Castrejón-Pita, A. A. 2019 An adaptive solver for viscoelastic incompressible two-phase problems applied to the study of the splashing of weakly viscoelastic droplets. J. Non-Newt. Fluid Mech. 264 (November 2018), 144–158.
  * Martínez Narváez et al. (2021) Martínez Narváez, C. D.V., Dinic, J., Lu, X., Wang, C., Rock, R., Sun, H. & Sharma, V. 2021 Rheology and Pinching Dynamics of Associative Polysaccharide Solutions. Macromolecules 54 (13), 6372–6388.
  * Mathues et al. (2018) Mathues, W., Formenti, S., McIlroy, C., Harlen, O. G. & Clasen, C. 2018 CaBER vs ROJER—Different time scales for the thinning of a weakly elastic jet. J. Rheol. 62 (5), 1135–1153.
  * McKinley (2005) McKinley, G. H 2005 Visco-elasto-capillary thinning and break-up of complex fluids. Rheol. Rev. 3 (05), 1–48.
  * Middleman (1965) Middleman, S. 1965 Stability of a viscoelastic jet. Chem. Eng. Sci. 20, 1037–1040.
  * Morrison & Harlen (2010) Morrison, N. F. & Harlen, O. G. 2010 Viscoelasticity in inkjet printing. Rheol. Acta 49 (6), 619–632.
  * Popinet (2003) Popinet, S. 2003 Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries. J. Comput. Phys. 190 (2), 572–600.
  * Popinet (2009) Popinet, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. J. Comput. Phys. 228 (16), 5838–5866.
  * Popinet (2018) Popinet, S. 2018 Numerical Models of Surface Tension. Annu. Rev. Fluid Mech. 50 (1), 49–75.
  * Rajesh et al. (2022) Rajesh, S., Thiévenaz, V. & Sauret, A. 2022 Transition to the viscoelastic regime in the thinning of polymer solutions. Soft Matter 18 (16), 3147–3156.
  * Rayleigh (1879) Rayleigh, L. 1879 On the Capillary Phenomena of Jets. Proc. R. Soc. London 29 (196-199), 71–97.
  * Renardy (2000) Renardy, M. 2000 Asymptotic structure of the stress field in flow past a cylinder at high Weissenberg number. J. Non-Newt. Fluid Mech. 90 (1), 13–23.
  * Renardy & Renardy (2004) Renardy, M. & Renardy, Y. 2004 Similarity solutions for breakup of jets of power law fluids. J. Non-Newt. Fluid Mech. 122 (1-3), 303–312.
  * Scharfman et al. (2016) Scharfman, B. E., Techet, A. H., Bush, J. W.M. & Bourouiba, L. 2016 Visualization of sneeze ejecta: steps of fluid fragmentation leading to respiratory droplets. Exp. Fluids 57 (2), 1–9.
  * Sharma et al. (2015) Sharma, V., Haward, S. J., Serdy, J., Keshavarz, B., Soderlund, A., Threlfall-Holmes, P. & McKinley, G. H. 2015 The rheology of aqueous solutions of ethyl hydroxy-ethyl cellulose (EHEC) and its hydrophobically modified analogue (hmEHEC): Extensional flow response in capillary break-up, jetting (ROJER) and in a cross-slot extensional rheometer. Soft Matter 11 (16), 3251–3270.
  * Snoeijer et al. (2020) Snoeijer, J. H., Pandey, A., Herrada, M. A. & Eggers, J. 2020 The relationship between viscoelasticity and elasticity. Proc. R. Soc. A Math. Phys. Eng. Sci. 476 (2243), 20200419\.
  * Tirtaatmadja et al. (2006) Tirtaatmadja, V., McKinley, H. G. & Cooper-White, J. J. 2006 Drop formation and breakup of low viscosity elastic fluids: Effects of molecular weight and concentration. Phys. Fluids 18 (4), 1–18.
  * Tomé et al. (2007) Tomé, M. F., Grossi, L., Castelo, A., Cuminato, J. A., McKee, S. & Walters, K. 2007 Die-swell, splashing drop and a numerical technique for solving the Oldroyd B model for axisymmetric free surface flows. J. Non-Newt. Fluid Mech. 141 (2-3), 148–166.
  * Turkoz et al. (2018) Turkoz, E., Lopez-Herrera, J. M., Eggers, J., Arnold, C. B. & Deike, L. 2018 Axisymmetric simulation of viscoelastic filament thinning with the Oldroyd-B model. J. Fluid Mech. 851, 1–13.
  * Turkoz et al. (2021) Turkoz, E., Stone, H. A., Arnold, C. B. & Deike, L. 2021 Simulation of impulsively induced viscoelastic jets using the Oldroyd-B model. J. Fluid Mech. 911, 1–29.
  * Villermaux (2007) Villermaux, E. 2007 Fragmentation. Annu. Rev. Fluid Mech. 39, 419–446.
  * Wagner et al. (2005) Wagner, C., Amarouchene, Y., Bonn, D. & Eggers, J. 2005 Droplet detachment and satellite bead formation in viscoelastic fluids. Phys. Rev. Lett. 95 (16), 164504.
  * Wagner et al. (2015) Wagner, C., Bourouiba, L. & McKinley, G. H. 2015 An analytic solution for capillary thinning and breakup of FENE-P fluids. J. Non-Newt. Fluid Mech. 218, 53–61.
  * Womersley (1955) Womersley, J. R. 1955 Method for the calculation of velocity rate of flow and viscous drag in arteries when the pressure gradient is known. J. Physiol. 127 (3), 553–563.
  * Xu et al. (2021) Xu, M., Li, X., Riseman, A. & Frostad, J. M. 2021 Quantifying the effect of extensional rheology on the retention of agricultural sprays. Phys. Fluids 33 (3).
  * Yapici et al. (2009) Yapici, K., Karasozen, B. & Uludag, Y. 2009 Finite volume simulation of viscoelastic laminar flow in a lid-driven cavity. J. Non-Newt. Fluid Mech. 164 (1-3), 51–65.
  * Yesilata et al. (2006) Yesilata, B., Clasen, C. & McKinley, G. H. 2006 Nonlinear shear and extensional flow dynamics of wormlike surfactant solutions. J. Non-Newt. Fluid Mech. 133 (2-3), 73–90.

(a)

(b)

Figure 15: Mesh convergence study showing the evolution of the minimum neck
radius $R_{min}(t)$ for a representative Lagrangian element with three
distinct levels-of-refinement $LVL=13,14$, and 15 for the parameters shown in
Table 2; results are presented through the temporal evolution of (a) the
minimum radius of the jet $R_{min}$, and (b), the local dimensionless strain
rate $Wi=\tau\dot{\epsilon}_{min}(t)$ to evaluate the resolution of the point
‘P3’ where the local dimensionless strain rate $Wi(t)$ passes through a
minimum.

##  Appendix A Mesh convergence study

It is essential to confirm that the current mesh resolution is adequate to
capture all the dynamics of interest, in particular, the point ‘P3’ labelled
in Figure 9(b) or equivalently the value of $Wi_{EC}$ which determines the
local rate of filament thinning in the elasto-capillary regime. To achieve the
required level of resolution we compute and present in Figure 15 the evolution
in $R_{min}(t)$ and $Wi(t)$ for a lower and a higher level of refinement,
$LVL=13$ and $LVL=15$, respectively. We note here that each increase of $LVL$
corresponds to a decrease in the minimum cell size by a factor of 2. For
example, $LVL=14$ and $LVL=15$ correspond to $\Delta x=0.006$ and $\Delta
x=0.003$, respectively. Figure 15(a) shows that there is a significant
influence of the grid resolution on the nonlinear elasto-capillary thinning
regime (region ‘IV’). Specifically, the local rate of thinning in the thread
radius for $LVL=13$ is excessively rapid and pinch-off is approached too
rapidly. The change in the slope is large as we move to a higher level of
resolution ($LVL=14$). However, increasing further to $LVL=15$ does not seem
to change the observed dynamics significantly, particularly regarding the
local exponential decrease of the minimum radius in the elasto-capillary
regime. An interesting difference between $LVL=14$ and $LVL=15$ is only
observed very close to the limit of the finite extensibility of the polymer
chains. A higher level of refinement delays the filament breakup and ensures a
slightly longer-lasting filament as the finite extensibility effects which
lead to the final breakup become significant later at $t-t_{p}\approx 8.4$
with $LVL=15$ compared to $t-t_{p}\approx 7.8$ with $LVL=14$. Additionally,
the enhanced resolution also results in a better-resolved terminal linear
thinning regime where finite extensibility effects dominate. This observation
can be confirmed by the evolution in the local Weissenberg number shown in
Figure 15(b). While $LVL=13$ refinement definitely leads to faster thinning
dynamics, it does not allow for meaningful quantitative analysis, the higher
resolution of the $LVL=14$ and $LVL=15$ simulations lead to good overlap for
the two minima‘min1’ and ‘min2’, in particular at point ‘P3’. Nonetheless, an
even closer inspection shows that the ‘P4’ point is not identical even at
these two levels of resolution. Therefore, $LVL=14$ refinement is sufficient
for computations of the exponential elasto-capillary regime, whereas $LVL=15$
refinement is required (but is also more computationally challenging) for
analysis of the terminal finite-extensibility regime that dominates the
ligament dynamics immediately before pinch-off.

##  Appendix B Linear stability analysis

Here, we provide details of the (temporal) linear stability analysis discussed
in Section 3.2 in connection with region ‘II’ that is identified in Figure 7.
Here we focus on the real part of the complex dispersion equation as the
analysis is restricted only to examining purely temporal instabilities of the
jet. Following the substitution of small amplitude perturbations in the
filament radius and the linearisation of the dimensionless mass and momentum
conservation equations in cylindrical coordinates, and the incorporation of
the kinematic and dynamic interfacial boundary conditions, the use of normal
mode analysis leads to the following dispersion relationship for an
axisymmetric Oldroyd-B jet, corresponding to a dilute polymer solution with
infinite chain extensibility, $L^{2}\rightarrow\infty$ Brenn et al. (2000):

$\displaystyle\Omega_{r}^{2}\frac{k}{2}\left[\frac{I_{0}\left(k\right)}{I_{1}\left(k\right)}+\frac{\rho_{g}}{\rho_{l}}\frac{K_{0}\left(k\right)}{K_{1}\left(k\right)}\right]+\Omega_{r}k^{2}Oh\frac{1+\beta
De\Omega_{r}}{1+De\Omega_{r}}$ (17)
$\displaystyle\quad\times\left[2k\frac{I_{0}\left(k\right)}{I_{1}\left(k\right)}\left(1+k^{2}\frac{Oh}{\Omega_{r}}\frac{1+\beta
De\Omega_{r}}{1+De\Omega_{r}}\right)\right.$
$\displaystyle\left.-1-2l\frac{I_{0}\left(l\right)}{I_{1}\left(l\right)}k^{2}\frac{Oh}{\Omega_{r}}\frac{1+\beta
De\Omega_{r}}{1+De\Omega_{r}}\right]$
$\displaystyle\quad=\frac{k^{2}}{2}\left(1-k^{2}\right)+C\frac{\rho_{g}}{\rho_{l}}k^{3}We\frac{K_{0}\left(k\right)}{K_{1}\left(k\right)}.$

Here, $\Omega_{r}>0$ ($\Omega_{r}<0$) indicates instability (stability),
$I_{n}$ and $K_{n}$ are the modified Bessel functions, $C$ is an empirical
correction factor to express the aerodynamic effects on the jet (here we
choose $C=0.175$ (Brenn et al., 2000; Keshavarz et al., 2015)), and $l$ is a
dimensionless modified wavenumber given by
$l^{2}=k^{2}+\left(1+De\left(\Omega+ik{U}_{0}\right)\right)/\left(Oh(1+\beta
De\left(\Omega+ikU_{0}\right))\right)$. The intrinsic Deborah number is
$De=\tau/t_{R}$, the Ohnesorge number is $Oh=\eta_{0}/\sqrt{\rho R_{0}\gamma}$
and the solvent viscosity ratio is $\beta=\eta_{s}/\eta_{0}$, with the
Rayleigh time $t_{R}=\sqrt{\rho R_{0}^{3}/\gamma}$ and the initial radius of
the jet $R_{0}$ being the characteristic time and length scales, respectively.
Eq. (17) can be solved numerically with a simple MATLAB-solver at a specific
value of $De$, $Oh$ and We for a range of wavenumbers with the dimensionless
growth rate $\Omega_{r}$ being the unknown. Typical results are presented in
Figure 7(b): In the inviscid limit the maximum growth rate $\Omega_{r}\approx
0.34$ is at $k\approx 0.693$. The first effect of viscosity is to stabilise
the jet partially and the most unstable growth rate reduces to
$\Omega_{r}\approx 0.23$ at $k\approx 0.6$. In contrast, linear viscoelastic
effects are observed to render the jet slightly more unstable than the
corresponding Newtonian viscous jet with the most unstable growth rate
increasing to $\Omega_{r}\approx 0.24$ at $k\approx 0.6$.