Transition to elasto-capillary thinning dynamics in viscoelastic jets

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, Gao & Bhat 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 atomization 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 optimization 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 the 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 (Turkoz et al. 2018; Snoeijer et al. 2020), highlighting the local importance of the polymeric extra 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, Entov & Rozhkov 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, Hinch & Li 2006; Bhat, Basaran & Pasquali 2008; Bhat et al. 2010; Morrison & Harlen 2010; Turkoz et al. 2018; Eggers, Herrada & Snoeijer 2020). The analytic studies have considered primarily infinitely extensible polymeric chains in a viscous solvent that can be described by the Oldroyd-B (Bird et al. 1987) constitutive equation (Li & Fontelos 2003; Clasen et al. 2006; Ardekani, Sharma & McKinley 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, Bourouiba & McKinley 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, Guan & Fu 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, Clasen & McKinley 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; Lauser, Rueter & Calabrese 2021; Martínez Narváez et al. 2021; Zinelis et al. 2024), as well as for more standard dripping experiments that focus on the transition to elasticity-dominated thinning (Amarouchene et al. 2001; Wagner et al. 2005; Rajesh, Thiévenaz & Sauret 2022).

Although the establishment of an elasto-capillary 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 that 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 flow rates due to the so-called ‘gobbling phenomenon’, with the conventional elasto-capillary balance expected in the CaBER instrument being established for larger flow rate 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 measuring accurately 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 flow rates. 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 incorporates the viscoelastic shearing flow upstream of the nozzle exit in addition to the subsequent capillarity-driven evolution of the jet and the formation of BOAS morphologies. Previous experimental work (Ghafourian et al. 1991; Mayer & Branam 2004; Lefebvre & McDonell 2017) has demonstrated that the growth of three-dimensional interfacial disturbances, and the transition from axisymmetric to three-dimensional structures is determined by the jet speed, gas-to-liquid density ratio and fluid elasticity (Liu & Liu 2008; Ruo, Chang & Chen 2008). Here, we study low-speed, moderately viscoelastic jets, for which the interfacial deformations are expected to remain axisymmetric (Bechtel, Cao & Forest 1992) (as confirmed with three-dimensional simulations provided in the supplementary material), by implementing the formulation of López-Herrera, Popinet & Castrejón-Pita (2019) whereby the azimuthal component of the extra polymeric stress tensor is also computed. We also follow these authors in solving a transport equation for the conformation tensor after the log-transformation (Fattal & Kupferman 2005) has been applied.

The rest of this paper is organized as follows. In § 2, the problem formulation and numerical procedure used to carry out the computations are outlined. In this section, we also validate our numerical predictions of the fully developed flow upstream of the nozzle exit plane against analytical solutions. A discussion of the results is provided in § 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 flow rate and the finite extensibility of the dissolved polymer. Finally, concluding remarks are provided in § 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 set-up 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$ that 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 § 2.2) 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

(2.1)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} =0, \end{gather}$$

(2.2)$$\begin{gather}\rho\left(\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}\right)={-}\boldsymbol{\nabla} p+\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma}+\gamma \kappa \boldsymbol{n} \delta_{s}, \end{gather}$$

where $t$, $\rho$, $\boldsymbol {u}$, $p$, $\boldsymbol {\sigma }$, $\gamma$, $\kappa$, $\boldsymbol {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 {\boldsymbol 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 } \boldsymbol {u}+(\boldsymbol {\nabla } \boldsymbol {u})^{\rm 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:

(2.3)\begin{equation} \boldsymbol{\sigma}_{p}=\frac{\eta_{p}} {\tau} \left(\frac{\boldsymbol{\mathsf{A}}}{1-\dfrac{\operatorname{tr}(\boldsymbol{\mathsf{A}})} {L^2}} -\boldsymbol{I} \right). \end{equation}

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 $\boldsymbol {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 equation

(2.4) \begin{equation} \frac{\partial \boldsymbol{\mathsf{A}}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\mathsf{A}} -(\boldsymbol{\nabla} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\mathsf{A}}+ \boldsymbol{\mathsf{A}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^{\rm T})={-}\frac{1}{\tau} \left(\frac{\boldsymbol{\mathsf{A}}}{1-\dfrac{\operatorname{tr}(\boldsymbol{\mathsf{A}})}{L^2}} -\boldsymbol{I} \right). \end{equation}

The above equations are rendered dimensionless by using the nozzle radius $R_0$, the Rayleigh velocity $U_R= \sqrt {\gamma / ( \rho _l R_0 )}$, 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 (2.1)–(2.4) leads to the following dimensionless equations:

(2.5)$$\begin{gather} {\tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}} =0, \end{gather}$$

(2.6)$$\begin{gather}\tilde{\rho} \left(\dfrac{\partial \tilde{\boldsymbol{u}}}{\partial \tilde{t}} + \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla} }\tilde{\boldsymbol{u}} \right) ={-} \tilde{\boldsymbol{\nabla}} \tilde{p} + Oh \left(\beta\, \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{\boldsymbol{\sigma}}_s + \frac{1-\beta}{De}\,\tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{\boldsymbol{\sigma}}_p \right) +\tilde{\kappa} \tilde{\delta} \boldsymbol{n}, \end{gather}$$

(2.7) $$\begin{gather} \frac{\partial \boldsymbol{\mathsf{A}}} {\partial \tilde{t}} + \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \tilde{\boldsymbol{\boldsymbol{\nabla}}} \boldsymbol{\mathsf{A}} -(\tilde{\boldsymbol{\boldsymbol{\nabla}}} \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\mathsf{A}}+ \boldsymbol{\mathsf{A}} \boldsymbol{\cdot} \tilde{\boldsymbol{\boldsymbol{\nabla}}} \tilde{\boldsymbol{u}}^{\rm T})={-} \frac{1}{De} \left(\frac{\boldsymbol{\mathsf{A}}}{1-\dfrac{\operatorname{tr}(\boldsymbol{\mathsf{A}})}{L^2}} -\boldsymbol{I}\right), \end{gather}$$

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} =-2\,{\rm D} (\log (R_{min}^{[\alpha ]})) /{\rm D}t$ 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,\ldots$ (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:

(2.8)\begin{equation} Wi ={-}2 \tau\,\frac{{\rm D}\log(R_{min}^{[\alpha]})}{{\rm D}t}.\end{equation}

This definition will be used to characterize 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, McKinley & Cooper-White 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 (2.3). This is known as the high-Weissenberg number problem (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 {\varPsi }$ of the conformation tensor $\boldsymbol {A}$, which can be computed via a diagonal transformation such that

(2.9)$$\begin{gather} \boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{R}} \boldsymbol{\varLambda} \boldsymbol{\mathsf{R}}^{\rm T}, \end{gather}$$

(2.10)$$\begin{gather}\boldsymbol{\varPsi} = \log \boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{R}}\,(\log\boldsymbol{\varLambda})\, \boldsymbol{\mathsf{R}}^T, \end{gather}$$

where $\boldsymbol {\varLambda }$ is the diagonal matrix containing the principal eigenvalues of $\boldsymbol {A}$. The numerical procedure for the solution of the appropriate transport equation for the conformation tensor $\boldsymbol {A}$ when the log-conformation transformation is applied, is provided by López-Herrera et al. (2019).

A volume-of-fluid 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

(2.11)\begin{equation} \frac{\partial c }{\partial \tilde{t}}+\tilde{\boldsymbol{u}} \boldsymbol{\cdot} \tilde{\boldsymbol{\boldsymbol{\nabla}}} c=0. \end{equation}

Using the volume-of-fluid formulation, which corresponds to a one-fluid approach to solving the two-phase flow, the density $\tilde{\rho}$ and viscosity $\tilde{\eta}$ are then respectively given by

(2.12)$$\begin{gather} \tilde{\rho} = c + (1-c)\,\frac{\rho_{g}}{\rho_l}, \end{gather}$$

(2.13)$$\begin{gather}\tilde{\eta} = c + (1-c)\,\frac{\eta_{g}}{\eta_l}, \end{gather}$$

where 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, 2021; López-Herrera et al. 2019) is used here to carry out the computations. A piecewise linear interface calculation 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 discretization 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 $100 R_0 \times 100 R_0$ as shown in figure 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. Downstream from the nozzle, a no-slip immersed boundary (as described in detail below) is applied for $r \geq 10 R_0$, which is far away from the jet region, which is located sufficiently far away from the jet region, for consistency with the immersed boundary method (Aniszewski et al. 2020) implemented to simulate the flow inside the nozzle as presented in details below. A pressure gradient is imposed at the left-hand boundary to drive the fluid through the nozzle, while atmospheric pressure is imposed at the right-hand 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}=A_{\theta \theta }=1$.

An adaptive mesh refinement 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 (Aniszewski et al. 2020) is employed to model the velocity field with no-slip and no-penetration conditions imposed at the solid walls. 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 a 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, Karasozen & Uludag 2009). As we show in figures 2(a,b); we obtain excellent agreement between the computed velocity and stress fields and the analytic results.

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 functions 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$.

Table 1. Parameter values utilized to validate the flow predictions within the nozzle using the immersed boundary method.

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 (2.14) with a dimensionless wavenumber $k=0.6$ (with tildes suppressed):

(2.14)\begin{equation} -{\boldsymbol{\nabla} } {p} \boldsymbol{\cdot} \boldsymbol{\hat{n}}_I = 8\, Ca\,(1+\epsilon_{p} \sin(\sqrt{We} \,k {t} )), \end{equation}

where $\hat{\boldsymbol{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 (2.6) through the relation $Ca = \sqrt {We} \, Oh$. The parameter values used in the simulations are provided in table 1. The value $k=0.6$ is selected according to the linear stability analysis presented in § 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)r\, {\rm d} r$ and axial velocity $\bar {u}_{x}=\int ^1_0u_x(r,x=0)r\,{\rm d} r$ 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)r\, {\rm d} r$ and $\bar {\sigma }_{p,rx}=\int ^1_0\sigma _{p,rx}(r,x=0)r\, {\rm d} r$, 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 at short times $t \leq 5$ after the flow is impulsively started. 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.

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 (2.14) with $\epsilon _p=0.4$ and $k=0.6$.

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$ ($0.02 > \Delta x_{minimum}>0.006$); 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 that develops upstream of the nozzle exit plane, as expected from figure 2(a), and discussed in § 2. This contour plot also shows that material elements at the centreline of the liquid jet leave the nozzle at the maximum velocity, and the velocity field rapidly rearranges (within one perturbation wavelength) so the fluid 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.

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, 2021). The dimensionless wavenumber $k$ is selected to be the value that corresponds to the most unstable mode of perturbation expected from linear stability analysis.

Figure 4. Simulation results at $t=47.8$ for the parameter values listed in table 2. (a) The volume fraction (left) and the mesh (right), 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 that 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 the polymeric stress components in the polymeric thread highlighting the highly localized distribution of the elastic stress and fore–aft asymmetry of the satellite drop that develops far downstream of the nozzle.

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 BOAS structures 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 adaptive mesh refinement 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 upper side of the thread (i.e. the side that 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 $\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.

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 {\rm \pi}n / (\sqrt {We} \, k)$, with $n=0,1,2,\ldots$). Mesh resolution represented by different $LVL$ values increases as time increases to ensure that the dynamics is 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 5b) 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 in both 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 detachment of the leading droplet, observed at early times, after which periodicity of the jet evolution is established. The spatio-temporal 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 5. (a) A sequence of jet profiles $R(x,t_i)$ highlighting the formation of BOAS 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 3a). (b) Space–time diagram or ‘kymograph’ showing $R_{min}(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.

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 characterized 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, Liu & Durst 2000) characterized by a sinusoidal response to the pressure gradient forcing set by (2.14).

Figure 6. Temporal evolution of the jet radius associated with four distinct regions. Regions I and II are 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 is strongly nonlinear, featuring the formation of the first satellite beads and BOAS 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.

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 enters into a fully nonlinear regime 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$) is 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 BOAS 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. regions I and II, before embarking on a detailed analysis of regions III and IV, in which the dynamics becomes increasingly nonlinear. In experimental visualization 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 us 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 BOAS structure develops fully (at $x \approx 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,\ldots$, 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. (a) 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 filled ($R_{min1}^{[\alpha ]}$, min1) and hollow ($R_{min2}^{[\alpha ]}$, min2) symbols of the corresponding colour for each neck region established behind each primary bead $[\alpha ]$. (b) 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.

3.2. Regions I and II: linear and weakly nonlinear evolution of jet profile

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

(3.1)\begin{equation} R_{min}(t) = 1 - \delta \exp(\varOmega t), \quad \text{for}\ t \geq 0,\end{equation}

where $\delta$ is the initial perturbation amplitude, and $\varOmega =\varOmega _r + {\rm i} \varOmega _i$ is the complex growth rate, which depends parametrically on $Oh$, $De$ and $k$; $\varOmega _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$ that 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, that is best characterized by the linear theory.

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 in the deviation of the jet radius from its base state for an imposed perturbation given by (3.1) 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 that 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 characterized 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 (B1) (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.

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 $\varOmega _r$ in (3.1). The growth rates obtained from the numerical simulation for $Oh=0.2$, $De=1$, $We=16$ and linear stability theory (Brenn et al. 2000) (see Appendix B for details) are $\varOmega _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, the finite extensibility of the polymeric chains is not expected to exert a strong influence on the flow at early times because the motion is driven by the action of capillarity and dominated by the linear interfacial disturbances. The finite extensibility of the polymer chains becomes critical as a thin thread is formed and undergoes severe thinning (Wagner et al. 2015), as will also be shown in § 3.4. We can therefore compare the numerically predicted growth rates with those obtained from linear stability theory for an Oldroyd-B fluid ($L^2 \rightarrow \infty$) 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 dash-dotted line), while the effects of viscosity 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 slightly 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 (i.e. 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 $\varOmega _r=0$ as obtained when $We=0$. In addition, the presence of fluid inertia is destabilizing, 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 $\varOmega _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,0.8$ (obtained using a similar procedure to that discussed in connection to figure 8a) shows good agreement. Specifically, each of 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$, $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.

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 they travel 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, characterized 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\,\%$ of $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 primary 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).

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 (EC) 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.

To determine the characteristic time scale that best describes the exponential thinning, the temporal evolution of the local dimensionless strain rate $Wi$, defined in § 2.1, is plotted in figure 9(b). Specifically, this plot highlights the existence of two different plateau values, $Wi=2/3$ and $Wi=1$, which correspond to the two distinct thinning rates during the elasto-capillary regime as determined by Keshavarz et al. (2015) and Mathues et al. (2018), respectively. In addition, figure 9(b) is characterized 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). Point 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 characterized by a local value of the Weissenberg number, which we will denote generically 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, BOAS 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 $Wi_{EC} \approx 1$ is established only for a narrow time interval, 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 figures 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 (3.1)) and the local strain rate $\dot {\epsilon }_{min} = -(2/R_{min})({\rm d}R_{min} / {\rm d}t)$ 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 centreline and the polymeric stress component at the jet centreline for two Lagrangian points corresponding to min1 and min2. Also shown 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). The centerline velocity associated with both min1 and min2, which correspond to the same Lagrangian element in the initial stages (regions I and 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 BOAS structure starts forming: the axial velocity associated with point min1 (min2) decreases (increases) until reaching a minimum (maximum). Shortly after reaching the 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.

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 simulation parameters provided in table 2. Also shown in (a) and (b) are profiles that depict the shape of the jet, coloured by contours indicating the magnitudes 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.

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 seen clearly. As expected, this results in the largest values of the axial component of the polymeric stress in the thin ligament (Clasen et al. 2006; Deblais et al. 2020; Eggers 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 BOAS 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 influence of the injection flow rate and the finite extensibility of the polymeric chains on the rate of thinning during the elasto-capillary regime and the associated pinch-off dynamics.

3.4. Effect of inertia and polymer chain finite extensibility

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 at $We \approx 2$; Clasen et al. 2009). It is seen that the jet length increases with $We$, and the thinning dynamics is 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, appears to be only 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.

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 the Weber number.

We also study the temporal evolution in the local dimensionless strain rate $Wi(t)=\tau \, \dot {\epsilon } _{min} (t)$ in figure 12, 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,\ldots$, we focus on only 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, and 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$), the transition from the characteristic points labelled P3 and P4 in figure 10 is more gradual compared to $We=16$ and $36$.

Figure 12. The effect of increasing the relative significance of fluid inertia (by increasing $We$) on the temporal evolution of the local dimensionless strain rate $Wi(t) = \tau \, \dot {\epsilon }_{min}(t)$ in the thinning neck of the fluid jet, for $We=8$, 16 and 36, with the rest of the parameters remaining unchanged from table 2.

In figure 13, we show how the dimensionless strain rate in a representative fluid neck varies with the extensibility parameter $L^2$ at 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_{EC}=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 9b), 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_{EC}=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.

Figure 13. The effect of increasing the extensibility parameter $L^2$ on the temporal evolution of the local dimensionless strain rate $Wi(t)= \tau \,\dot {\epsilon }_{min}(t)$, for $L^2=900, 2500, 10\,000$ and $We=16$, with the rest of the parameters remaining unchanged from table 2. Also shown for comparison is the Oldroyd-B limit corresponding to $L^2 \rightarrow \infty$.

In figure 14, we construct a flow map in $(We,L^2)$ 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 $(We,L^2)$ 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 that highlight the shape of the jet for low and high values of $We$ and $L^2$ are coloured by the magnitude of the axial component of the elastic stress. The green triangles correspond to the numerical simulations that have been performed here to capture how the magnitude of $Wi_{EC}$ varies in the $(We,L^2)$ space. 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 BOAS 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 characterized by a large value of the polymer finite extensibility parameter, an elasto-capillary balance with $Wi_{EC} =2/3$ (deep red colours) is established, similar to the dynamics realized in the CaBER device (Entov & Hinch 1997; Anna et al. 2001). However, when the axial momentum in the jet is high, the length to breakup is large, and the finite extensibility of the polymeric chains is small, the asymmetric force balance of Clasen et al. (2009) and Mathues et al. (2018) applies, resulting in a faster local stretching rate such that $Wi_{EC} \approx 1$ (yellow colour contours).

4. Conclusions

We have studied the thinning and breakup of an axisymmetric viscoelastic jet issuing from a nozzle using the FENE-P constitutive relation that 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 utilizes 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 (BOAS) 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 characterized by the formation of a distinct BOAS 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 $R_{min} \approx 0.006 R_0$ and jet lengths as large as $\ell _{jet,max} \approx 100 R_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 BOAS 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 Weber number and chain extensibility space, and calculate the variations in the local dimensionless extension rate $Wi_{EC} (We, L^2)$. This map helps us to identify regions of parameter space characterized 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$ to $Wi_{EC} =1$. 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 inkjet device (Morrison & Harlen 2010; Keshavarz et al. 2015; Mathues et al. 2018; Xu et al. 2021).