2.1. General two-dimensional formulation

In what follows, variables or parameters with tilde (˜) over them, are denoted as dimensional while their dimensionless counterparts are presented without tilde. We consider a 2-D, axisymmetric liquid bridge of an EVP material held captive between two parallel and coaxial disks. The initial separation of the two disks is $\tilde{L}_{0}$ and their common radius is $\tilde{h}_{0}$ , both of which are of the order of a few millimetres, at the most, making gravitational effects very small from the beginning of the process. The material is pinned to the disk’s perimeter, so that the initial radius of the fluid thread is also $\tilde{h}_{0}$ . At time $\tilde{t}=0^{+}$ , the upper plate starts to move with constant velocity in the positive $\tilde{z}$ -direction. As time proceeds, the capillary force and material stress components start to grow, and a neck is formed in the filament. We monitor the time evolution of the minimum radius of the filament as it shrinks to zero, $\tilde{h}_{\textit{min}}\rightarrow 0$ , along with the dominant velocity and stress components. The negligible gravity force allows us to assume a symmetry plane at $\tilde{z}=0$ , making it possible to simulate only the top half of the filament. Figure 1 shows a schematic of the axisymmetric EVP bridge being stretched with symmetry plane imposed at $\tilde{z}=0$ .

Figure 1. Schematic of the axisymmetric EVP filament stretching geometry.

The characteristic scales used for making the variables dimensionless are those proposed by Moschopoulos et al. (Reference Moschopoulos, Kouni, Psaraki, Dimakopoulos and Tsamopoulos2023) for the bulk dynamics of an EVP filament stretching. The radius of the disk $\tilde{h}_{0}$ and the capillary force $({\tilde{\sigma }}/{\tilde{h}_{0}})$ , scale all lengths and pressure or stress components, respectively. Balancing capillarity and viscous stresses results in the visco-capillary time, which is used as the characteristic time scale, $\tilde{t}_{vc}=({\tilde{k}\tilde{h}_{0}}/{\tilde{\sigma }})^{({1}/{n})}$ , where $\tilde{k}\, ({\textrm{Pa s}}^{n})$ is the consistency index, and $n$ is the shear-thinning exponent of the Saramito–Herschel–Bulkley (SHB) model. The velocity scale is ${\tilde{h}_{0}}/{\tilde{t}_{vc}}$ . Another form of stress scale, equivalent to the one given earlier, is ${\tilde{k}}/{{(\tilde{t}_{vc}})^{n}}$ . Introducing these characteristic scales in the governing equations results in the dimensionless numbers for the present analysis, which are given in table 1.

Table 1. Dimensionless numbers.

With these characteristic quantities, the dimensionless mass and momentum conservation equations become

(2.1) \begin{align}& Oh^{-2}\frac{D\boldsymbol{u}}{Dt} =-\boldsymbol{\nabla }p+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }_{\boldsymbol{s}}-Bo\, \boldsymbol{e}_{\boldsymbol{z}}, \end{align}

(2.2) \begin{align}&\qquad\qquad\qquad\qquad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} =0, \end{align}

where ${\textrm{D}}/{{\textrm{D}}t}$ is the substantial derivative, $\boldsymbol{u}$ the velocity vector, $\boldsymbol{\nabla}$ the gradient operator, $\boldsymbol{\tau }_{\boldsymbol{s}}$ the Newtonian solvent stress, $\boldsymbol{\tau }$ the polymeric extra stress tensor, p the pressure and $\boldsymbol{e}_{\boldsymbol{z}}$ the unit vector in the positive z-direction. In this study, the Newtonian solvent contribution is taken to be zero, hence $\boldsymbol{\tau }_{\boldsymbol{s}}=\textbf{0}$ . Furthermore, in this study we will neglect inertial effects, setting $Oh^{-2}=0$ . The polymeric extra stress tensor is described by the SHB (Saramito Reference Saramito2009) EVP model for the rheological response of the material

(2.3) \begin{equation}Ec\, \overset{\boldsymbol{\nabla }}{\boldsymbol{\tau }}+\max \!\left(0,\left| \boldsymbol{\tau}_{\boldsymbol{d}}\right| -Y_{s}\right)^{\frac{1}{n}}\frac{\boldsymbol{\tau }}{\left| \boldsymbol{\tau}_{\boldsymbol{d}}\right| }=\dot{\boldsymbol{\gamma }}.\end{equation}

The first term in this model is the elastic contribution and the second one is the viscoplastic contribution to the total rate of deformation tensor, $\dot{\boldsymbol{\gamma }}=\boldsymbol{\nabla }\boldsymbol{u}+(\boldsymbol{\nabla }\boldsymbol{u})^{T}$ . In the above model, $\overset{\boldsymbol{\nabla }}{\boldsymbol{\tau }}$ is the upper-convected derivative of extra stress tensor, which is defined as

(2.4) \begin{equation}\overset{\boldsymbol{\nabla }}{\boldsymbol{\tau }}=\frac{\partial \boldsymbol{\tau }}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\tau }-\left(\boldsymbol{\nabla }\boldsymbol{u}\right)^{T}\boldsymbol{\cdot }\boldsymbol{\tau }-\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}. \end{equation}

Also, $\boldsymbol{\tau}_{\boldsymbol{d}}=\boldsymbol{\tau}-(tr(\boldsymbol{\tau })/tr(\boldsymbol{I}))$ is the deviatoric part of extra stress tensor with magnitude $| \boldsymbol{\tau} _{\boldsymbol{d}}| =\sqrt{({1}/{2})\boldsymbol{\tau }_{\boldsymbol{d}}\colon \boldsymbol{\tau }_{\boldsymbol{d}}}$ , and $(\boldsymbol{\nabla }\boldsymbol{u})^{T}$ indicates the transpose of the velocity gradient. The max term determines the transition from a solid-like to a fluid-like region by incorporating the von Mises criterion. If $| \boldsymbol{\tau}_{\boldsymbol{d}}|$ exceeds $Y_{s}$ , the material yields and deforms like a viscoelastic fluid with shear thinning, otherwise it behaves like a neo-Hookean solid. The free surface of the filament, $S_{f}$ , is in contact with air, which is assumed to have negligible viscosity and density. In our formulation, nodes on $S_{f}$ are treated in a Lagrangian sense, hence the velocity of the mesh and the velocity of the fluid particle on $S_{f}$ in the direction normal to $S_{f}$ must be equal. This is the so-called kinematic boundary condition and serves as a boundary condition on the free surface. Furthermore, the total stress from the fluid side acting on $S_{f}$ must balance capillarity, given that air is dynamically inert. This is the so-called traction boundary condition. These two boundary conditions are

(2.5) \begin{align} \boldsymbol{n}\boldsymbol{\cdot }\left(\boldsymbol{u}-\boldsymbol{u}_{m}\right) = 0,\quad & {\textit{on}}\, S_{f}, \end{align}

(2.6) \begin{align} \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{T} =\left(2\mathcal{H}\right)\boldsymbol{n},\quad & {\textit{on}}\, S_{f}, \end{align}

where $\boldsymbol{u}_{m}$ is the velocity of the mesh nodes and $\boldsymbol{n}$ the outward unit vector normal to $S_{f}$ ; T the total stress tensor, defined as $\boldsymbol{T}=-p\boldsymbol{I}+\boldsymbol{\tau }$ , with $\boldsymbol{I}$ the unit tensor; $2\mathcal{H}=-\boldsymbol{\nabla} _{s}\boldsymbol{\cdot }\boldsymbol{n}$ is twice the mean curvature of $S_{f}$ and lastly $\boldsymbol{\nabla} _{s}=(\boldsymbol{I}-\boldsymbol{nn})\boldsymbol{\cdot }\boldsymbol{\nabla} \text{is the surface gradient operator}$ (Deen Reference Deen2012). We have assumed axial symmetry with axis at $r$ = 0, or on $S_{A}$ and planar symmetry at $z=0$ , or on $S_{p}$ . On these boundaries the no-shear and no-penetration boundary conditions are applied

(2.7) \begin{align} \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{T}\boldsymbol{\cdot }\boldsymbol{t} = 0, &\quad {\textrm{on}}\,S_{A}\, {\textrm{and}}\, S_{P}, \end{align}

(2.8) \begin{align} \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{u} = 0, &\quad {\textrm{on}}\,S_{A}\, {\textrm{and}}\, S_{P}, \end{align}

where $\boldsymbol{t}$ is the unit vector tangential to either $S_{A}$ or $S_{p}$ . The fluid in contact with the two solid disks, initially at rest, follows the no-slip and no-penetration boundary conditions and when t > 0, the upper solid disk starts to move suddenly with dimensionless velocity U

(2.9) \begin{align}& \begin{cases} t=0, \quad {\boldsymbol{u}}=0\\ t \gt 0, \quad {\boldsymbol{u}}=U{\boldsymbol{e}}_{\boldsymbol{z}} \end{cases} \text{at}\quad z=L, \end{align}

(2.10) \begin{align}&\qquad\quad \boldsymbol{u}=\textbf{0}, \quad \text{at} \quad z=0,\end{align}

where $L=L_{0}+Ut$ is the disk separation at time t. The contact point on the intersection between the free surface $S_{f}$ and the solid disks is pinned

(2.11) \begin{equation}h_{0}=1, \quad \text{at}\quad z=L\, \quad {\textrm{and}}\, \quad z=0.\end{equation}

Due to the hyperbolic nature of the SHB model (and similarly for most other constitutive laws for polymeric liquids), the boundary conditions for the extra stress tensor are needed only at the inflow (Van Der Zanden & Hulsen Reference Van Der Zanden and Hulsen1988). In this problem, no inflow conditions exist, and no stress boundary conditions are required, except for $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{T}\boldsymbol{\cdot }\boldsymbol{t}=0$ at the axis of symmetry and the symmetry plane along with the force balance at the free interface.

2.2. Slender thread equations

Assuming that to highest order the variation of all flow variables in the radial direction, $r$ , is negligible compared with that in the axial direction, $z$ (except for $u_{r},\, {\textrm{and}}\, \tau _{rz}$ , which vary linearly with $r$ ), the interface $S_{f}$ can be described as a function of $z$ alone, $S_{f}, r=h(z,t)$ , disallowing any overturning of the free surface. Then the previously presented 2-Dl equations can be reduced to a much simpler, 1-D equation set. These are the so-called slender filament equations. The reduced mass and momentum equations, still neglecting solvent viscosity, are (Bechtel et al. Reference Bechtel, Cao and Forest1992; Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006)

(2.12) \begin{align}&\qquad\qquad\qquad\quad \frac{\partial h^{2}}{\partial t}+\frac{\partial }{\partial z}\big(u_{z}h^{2}\big)=0, \end{align}

(2.13) \begin{align}& Oh^{-2}\left(\frac{\partial u_{z}}{\partial t}+u_{z}\frac{\partial u_{z}}{\partial z}\right)=\frac{\partial }{\partial z}(2\mathcal{H})+\frac{1}{h^{2}}\frac{\partial }{\partial z}\big(h^{2}(\tau _{\textit{zz}}-\tau _{\textit{rr}})\big), \end{align}

(2.14) \begin{align}&\qquad\qquad\qquad 2\mathcal{H}=\frac{h_{\textit{zz}}}{\left(1+h_{z}^{2}\right)^{\frac{3}{2}}}-\frac{1}{h\left(1+h_{z}^{2}\right)^{\frac{1}{2}}}, \end{align}

where $u_{z}(z,t)$ is the axial velocity, $2\mathcal{H}$ is twice the mean curvature of the free surface, now explicitly written in the cylindrical coordinate system and the subscript $z$ next to $h$ , as in $h_{z}$ , indicates partial differentiation; in this case, $({\partial h}/{\partial z}).$ Please note that in § 2.2 all variables have retained their usual symbols but depend only in the axial direction $z$ and time, $t$ . Equation (2.13) is the reduced form of the dominant momentum balance in the axial direction, in which pressure has been eliminated in favour of the curvature term via the normal force balance and the next-order term of the axial velocity has been eliminated via the tangential force balance. As stated in the previous section, we will analyse the slender filament equations ignoring inertia. Although the computational effort to solve the 1-D equations is considerably smaller than the effort to solve their 2-D counterparts, the former equations are not accurate enough to capture the flow field under certain conditions (Suryo & Basaran Reference Suryo and Basaran2006; Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018). In this study, we avoid conducting 1-D simulations; we will use the slender filament equations only to derive the asymptotic decay of the radius and growth of the stress in the pinch-off regime.

4.1. Asymptotic analysis near filament pinching at the symmetry plane via the slender thread equations following the SHB model

Since we are seeking a self-similar solution in the vicinity of the singular point, where $h_{\textit{min}}\rightarrow 0$ , we introduce two new local variables, $Z$ and $T$ , which are defined as

(4.1) \begin{align} Z&=z-z_{0}, \end{align}

(4.2) \begin{align} T&=t_{0}-t. \end{align}

In the above expressions, $Z$ denotes the axial distance from the singular point (the location of $h_{\textit{min}})$ , the axial coordinate of which is $z_{0}$ . Having set $Oh^{-2}=0$ , to neglect inertia in the momentum equation, we find $z_{0}=0$ in all our simulations and $z$ is probed where $h=1.01h_{\textit{min}}$ . Here, $T$ is the time until the breakup time, which is denoted as $t_{0}$ . Preliminary 2-D results indicate that the time evolution of the minimum radius follows the pattern shown in figure 2.

Figure 2. The evolution of minimum radius vs. time, for parameter values $n=0.45,\, Ec=0.14,\, Y_{s}=1.25,\, U=0.39$ .

This figure clearly demonstrates that the breakup time of the filament is finite, a trend observed in all our simulations for values of $n$ less than unity. This should have been expected, because the SHB model predicts a finite extensible material, which is known to have a finite breakup time as opposed to the Oldroyd-B fluid, which approaches breakup exponentially (Renardy Reference Renardy1995; McKinley Reference McKinley2005). Hence, we introduce self-similar forms with power-law functions for the axial position, axial velocity and axial component of extra stress tensor

(4.3) \begin{align} h\left(Z,T\right)&=T^{{\alpha _{1}}}\, \phi \left(\xi \right), \end{align}

(4.4) \begin{align} u_{z}\left(Z,T\right)&=T^{{\alpha _{2}}}\, \psi \left(\xi \right), \end{align}

(4.5) \begin{align} \tau _{\textit{zz}}\left(Z,T\right)&=T^{{\alpha _{3}}}{\Omega} \left(\xi \right), \end{align}

(4.6) \begin{align} \xi &=\frac{Z}{T^{\delta }}. \\[6pt] \nonumber \end{align}

Please note that in § 4 all variables have retained their previous symbols but now are functions of $Z,T$ . Moreover, $\xi$ denotes the similarity variable and the exponent $\delta$ controls the time dependence of the extent of the similarity region. When the filament radius becomes extremely small, both $u_{z}(Z,T)$ and $\tau _{\textit{zz}}(Z,T)$ in the vicinity of the singular point grow increasingly fast. Under these conditions, the nonlinear terms in the EVP constitutive equation become dominant. From now on, we will refer to this regime as the highly nonlinear regime. Next, we substitute the self-similar forms of $h(Z,T)$ and $u_{z}(Z,T)$ , in the continuity equation, (2.12), and enforce time invariance to obtain

(4.7) \begin{equation}\delta -\alpha _{2}=1,\quad \phi '=\phi \frac{\alpha _{1}-\dfrac{\psi '}{2}}{\psi +\delta \xi }.\end{equation}

Given that the axial momentum equation is the dominant momentum balance and $\tau _{\textit{zz}}$ the dominant stress component, we analyse the $z$ -component of the constitutive equation, which according to the SHB model is

(4.8) \begin{equation}\frac{\partial \tau _{\textit{zz}}}{\partial t}+u_{z}\frac{\partial \tau _{\textit{zz}}}{\partial z}=2\frac{\partial u_{z}}{\partial z}\tau _{\textit{zz}}+\frac{2}{Ec}\frac{\partial u_{z}}{\partial z}-\frac{\max \!\left(0,\, \left| \boldsymbol{\tau} _{\boldsymbol{d}}\right| -Y_{s}\right)^{\frac{1}{n}}}{Ec}\frac{{\tau} _{\textit{zz}}}{\left| \boldsymbol{\tau} _{\boldsymbol{d}}\right| }.\end{equation}

According to its definition, $| \boldsymbol{\tau}_{\boldsymbol{d}}|$ is

(4.9) \begin{equation}\left| \boldsymbol{\tau} _{\boldsymbol{d}}\right| =\sqrt{\frac{1}{2}\left(\left(\tau _{\textit{rr}}-\frac{tr\left(\boldsymbol{\tau} \right)}{3}\right)^{2}+2\tau _{rz}^{2}+\left(\tau _{\textit{zz}}-\frac{tr\left(\boldsymbol{\tau} \right)}{3}\right)^{2}+\left(\tau _{\theta \theta }-\frac{tr\left(\boldsymbol{\tau} \right)}{3}\right)^{2}\right)}.\end{equation}

In regions of high nonlinearity $tr(\boldsymbol{\tau })\approx \tau _{\textit{zz}}$ , yielding

(4.10) \begin{equation}\left| \boldsymbol{\tau} _{\boldsymbol{d}}\right| \approx \sqrt{\left(\frac{1}{2}\right)\left[\left(-\frac{\tau _{\textit{zz}}}{3}\right)^{2}+\left(\frac{2\tau _{\textit{zz}}}{3}\right)^{2}+\left(-\frac{\tau _{\textit{zz}}}{3}\right)^{2}\right]}=\sqrt{\frac{1}{3}\tau _{\textit{zz}}^{2}}=\frac{1}{\sqrt{3}}\left| \tau _{\textit{zz}}\right|.\end{equation}

Moreover, in the same region $| \tau _{\textit{zz}}| \gg Y_{s}$ , while in the upper half of the thread that we examine $\tau _{\textit{zz}}\geq 0$ , simplifying the max term in the EVP constitutive equation to $({1}/{\sqrt{3}}\tau _{\textit{zz}})^{{1}/{n}}$ . Therefore, the YS does not enter the dominant balance in the constitutive law and the entire thread is yielded in all examined cases. By assuming that the normal $z$ -component of the stress tensor in the EVP thread is constant, we neglect the convective part of the substantial derivative in (4.8). The rest of the terms in the EVP constitutive equation are approximated following ideas in Renardy (Reference Renardy1995) and Fontelos & Li (Reference Fontelos and Li2004), resulting in

(4.11) \begin{equation}\frac{\partial \tau _{\textit{zz}}}{\partial t}=2\frac{\partial u_{z}}{\partial z}\tau _{\textit{zz}}-\frac{\left(\frac{1}{\sqrt{3}}\tau _{\textit{zz}}\right)^{\frac{1}{n}}}{Ec}\sqrt{3}.\end{equation}

Substituting the assumed forms for the similarity variables and enforcing time invariance yields

(4.12) \begin{equation}\alpha _{3}=\frac{n}{n-1},\,\, -\alpha _{3}{\Omega} +\delta \xi {\Omega}^{\prime} =2\phi^{\prime}{\Omega} -\frac{1}{Ec\left(\sqrt{3}\right)^{n-1}}{\Omega} ^{\frac{1}{n}}.\end{equation}

Considering inertialess physics, the slender momentum equation, (2.13), becomes

(4.13) \begin{equation}\frac{\partial }{\partial z}\big(h+h^{2}(\tau _{\textit{zz}}-\tau _{\textit{rr}})\big)=0.\end{equation}

For this balance to be consistent, $\tau _{\textit{zz}}$ must grow like 1/h so we can conclude that

(4.14) \begin{equation}\alpha _{1}=-\alpha _{3}=\frac{n}{1-n}.\end{equation}

Turning our attention again to the mass conservation equation, we obtain

(4.15) \begin{equation} \phi^{\prime}=\phi \frac{\alpha _{1}-\dfrac{\psi^{\prime}}{2}}{\psi +\delta \xi }.\end{equation}

The denominator will be zero where $\psi =-\delta \xi$ , for $\xi =\xi _{0}$ . Therefore, as in the Newtonian (Doshi et al. Reference Doshi, Suryo, Yildirim, McKinley and Basaran2003; Eggers Reference Eggers1993) or the power-law material (Doshi et al. Reference Doshi, Suryo, Yildirim, McKinley and Basaran2003), a smooth solution exists only if, at $\xi =\xi _{0}$ , the following hold:

(4.16) \begin{align} \psi \left(\xi _{0}\right)&=-\delta \xi _{0}, \end{align}

(4.17) \begin{align} \psi^{\prime} \left(\xi \right)&=2\alpha _{1}. \end{align}

Hence, a Taylor series expansion in $\xi$ about $\xi =\xi _{0}$ as in Doshi et al. (Reference Doshi, Suryo, Yildirim, McKinley and Basaran2003) results in

(4.18) \begin{align} \phi \left(\xi -\xi _{0}\right)&=\phi _{0}+\phi _{2}\left(\xi -\xi _{0}\right)^{2}+\phi _{4}\left(\xi -\xi _{0}\right)^{4}+\ldots, \end{align}

(4.19) \begin{align} \psi \left(\xi -\xi _{0}\right)&=-\delta \xi _{0}+2\alpha _{1}\left(\xi -\xi _{0}\right)+\psi _{3}\left(\xi -\xi _{0}\right)^{3}+\ldots, \end{align}

where $\phi _{0},\, \phi _{2},\, \ldots $ and $\psi _{3},\, \ldots$ are the coefficients of the expansion. In all our simulations, the value of the axial velocity is zero at the plane of symmetry, which arises at $\xi _{0}=0$ for inertialess flow physics. Therefore, $\phi$ and $\psi$ are even and odd functions, respectively. The minimum value of $\phi$ is $\phi _{0}$ . Hence, the minimum radius takes the following power-law form:

(4.20) \begin{equation}h_{\textit{min}}\left(t\right)=\phi _{0}\, T^{{\alpha _{1}}},\end{equation}

where $\phi _{0}$ depends on physical properties and can be determined numerically. In this study, because of the dual nature of the material, and particularly the max term in the constitutive model, the flow is more complicated, and we do not determine analytically the value of $\phi _{0}$ . This could be pursued in future studies. The main focus of this study is the breakup asymptotic forms and the existence of self-similarity. In a nutshell, the results for the scalings in the capillary breakup regime are given by (4.7a ) and (4.14).

The asymptotes derived above are similar with those derived by Renardy, (Reference Renardy2002a ) for a generalised form of the viscoelastic Phan–Thien–Tanner (G-PTT) model. The G-PTT fluid described in Renardy, (Reference Renardy2002a ) in dimensional form can be written as follows:

(4.21) \begin{equation}\tilde{\lambda }\overset{\boldsymbol{\nabla }}{\tilde{\boldsymbol{\tau }}}+\tilde{\boldsymbol{\tau }}+\frac{\nu \tilde{\lambda }^{a-1}}{\tilde{\eta }_{p}^{a-1}}\left(tr\left(\tilde{\boldsymbol{\tau }}\right)\right)^{a-1}\tilde{\boldsymbol{\tau }}=\tilde{\eta }_{p}\widetilde{\dot{\boldsymbol{\gamma }}},\end{equation}

where $a\gt 1$ is the exponent of the G-PTT model, $\tilde{\lambda }$ is the relaxation time, $\tilde{\eta }_{p}$ is the polymeric viscosity and $\nu$ is a positive dimensionless number, playing the same role as the mobility factor in the Giesekus model. The $z$ -component of G-PTT model where $\tau _{\textit{zz}}$ is the dominant stress component is

(4.22) \begin{equation}\frac{\partial \tilde{\tau }_{\textit{zz}}}{\partial \tilde{t}}=2\frac{\partial \widetilde{u_{z}}}{\partial \tilde{z}}\tilde{\tau }_{\textit{zz}}-\frac{\nu \tilde{\lambda }^{a-2}}{\tilde{\eta }_{p}^{a-1}}\tilde{\tau }_{\textit{zz}}^{a}.\end{equation}

Equation (4.22) reduces to (4.11) if the exponent $a$ in (4.22) is replaced by $({1}/{n})$ . Hence, it is expected that the G-PTT model will have the same asymptotes as the SHB model in the breakup regime. According to Renardy, (Reference Renardy2002a ), the asymptotes of the G-PTT model are

(4.23) \begin{align} \alpha _{1}&=\frac{1}{a-1}, \quad \alpha _{2}=-1, \quad \alpha _{3}=\frac{1}{1-a}, \!\!\qquad\text{if }1\lt a\lt \frac{7}{3}, \end{align}

(4.24) \begin{align} \alpha _{1}&=\frac{1}{4a-8}, \quad \alpha _{2}=-1, \quad \alpha _{3}=\frac{1}{1-a}, \quad \text{ if }\, \frac{7}{3}\lt a\lt 3. \end{align}

For $({7}/{3})\lt a\lt 3\, (0.33\lt n\lt 0.43),$ surface tension has no influence on breakup, and pinching enters a new regime, the so-called elastic breakup regime. Our numerical results are consistent with this finding, as we will demonstrate in § 5.5.2. Integration of (4.13) yields

(4.25) \begin{equation}h+h^{2}\tau _{\textit{zz}}=F\left(t\right),\end{equation}

where $F(t)$ is the total force acting on the filament cross-section. In the elastic breakup regime, the first term on the left-hand side is associated with capillarity and is subdominant. Substituting $\alpha _{1}$ and $\alpha _{3}$ , as given in (4.24), into the asymptotic forms of (4.3) and (4.5), one obtains

(4.26) \begin{equation}F\left(t\right)\sim h^{2}\tau _{\textit{zz}}\sim \, \left(t_{0}-t\right)^{{2\alpha _{1}}+{\alpha _{3}}}\sim \left(t_{0}-t\right)^{\frac{2n-6n^{2}}{\left(n-1\right)\left(4-8n\right)}}.\end{equation}

For the force to be finite, requires that $({2n-6n^{2}})/{(n-1)(4-8n)}\gt 0$ which yields $n\gt {1}/{3}$ . Indeed, our code did not converge for values $n\lt 0.33$ . Furthermore, Renardy suggests that $\delta =0,\, \alpha _{2}=-1$ , but our 2-D numerical results slightly deviate from this value. The reason for this difference is probably caused by the fact that Renardy has neglected any radial dependence of $u_{z}, p,\, \tau _{\textit{zz}},\, \tau _{rr\, }\, {\textrm{and}}\, \tau _{\theta \theta }$ and completely neglected $u_{r}, \tau _{rz}$ , which is not entirely correct, as will be seen in § 5. On the other hand, our solution of the 2-D equations may include numerical error. We calculate the following scaling exponents for axial velocity and axial position numerically:

(4.27) \begin{equation}\alpha _{2}=\frac{n}{n-1}\left(\lt 0\right),\, \delta =\frac{n}{n-1}+1(\lt 0), \quad \textrm{if}\quad n\gt 0.5,\end{equation}

where and

(4.28) \begin{equation}\alpha _{2}\rightarrow -1, \delta \rightarrow 0,\quad \textrm{if}\quad n\leq 0.5\, \quad {\textrm{and}} \quad n\rightarrow 0.4.\end{equation}

In late stages of breakup, since $\tau _{\textit{zz}}$ is dominant over $\tau _{\textit{rr}}$ , the asymptotic behaviour of the transient extensional viscosity $\eta _{ext}=(\tau _{\textit{zz}}-\tau _{\textit{rr}})/\dot{\epsilon }_{\textit{local}}$ yields

(4.29) \begin{equation}\eta _{ext}\sim \left(t_{0}-t\right)^{\frac{2n-1}{n-1}}.\end{equation}

Clearly, for $n\gt 0.5$ , $\eta _{ext}$ diverges to infinity, while for $n\lt 0.5$ , $\eta _{ext}$ tends to 0, indicating transient extension-rate-thickening and -thinning behaviours, respectively. It is noteworthy that the exponents in EVP fluids are different from those in power-law or viscoplastic fluids. Similar deviations have been observed in experiments with shear-thinning viscoelastic materials where scalings do not follow expected asymptotes of a power-law material (Moon et al. Reference Moon, Yamani, McKinley and Lee2024). In all our simulations, the scaling exponent for the axial position, denoted as δ, is found to be negative. We are probing the axial position at $z$ where $h=1.01h_{\textit{min}}$ . The negative sign of δ implies that the neck length increases as breakup is approached, instead of decreasing, as is usually the case with $\delta \gt 0$ . Moreover, the magnitude of the scaling exponent of the minimum radius, |α₁|, consistently exceeds |δ| in our simulations. Hence, the thinning of the thread occurs faster than the growth of the length of the neck. Consequently, although breakup still occurs at $z=0$ , the axial extent of fast thinning increases. This has been called ‘breakup at a finite axial length’ by Renardy, (Reference Renardy2002a ) who found that $\delta =0$ , analytically for the G-PTT fluid. This implies that the neck length remains constant. Consequently, the inertialess breakup of EVP fluids occurs over a finite length, whereas in power-law or viscoplastic fluids, it takes place at a single point. These ideas will be discussed and supported further in § 5.5.4.

4.2. Ideal uniaxial extension of an EVP material under highly nonlinearity conditions

To further validate the derived scaling exponents for $h$ and $\tau _{\textit{zz}}$ , we solve numerically an ideal uniaxial extension flow with the EVP constitutive equation. This simplification reduces the slender jet equations to those for a perfectly cylindrical liquid jet under extension. Thus, the strain rate becomes

(4.30) \begin{equation}\dot{\epsilon }=\frac{\partial u_{z}}{\partial z}=-\frac{2}{h}\frac{{\textrm{d}}h}{{\textrm{d}}t}.\end{equation}

Assuming the power-law asymptote for the radius, the local strain rate reduces to

(4.31) \begin{equation}h=\phi _{0}\left(t_{0}-t\right)^{{\alpha _{1}}}\, \Longrightarrow \, \dot{\epsilon }=\frac{2\alpha _{1}}{t_{0}-t}\approx \frac{2}{T}.\end{equation}

We note that assuming $({2\alpha _{1}})/({t_{0}-t})\approx {2}/{T}$ will not affect the high nonlinearity trend. In uniaxial extension

(4.32) \begin{equation}u_{r}=-0.5\dot{\epsilon }r, \quad u_{\theta }=0, \quad u_{z}=\dot{\epsilon }z.\end{equation}

In this ideal flow, the velocity gradient $\boldsymbol{\nabla }\boldsymbol{u}$ and $\boldsymbol{\dot{\gamma}}$ are

(4.33) \begin{equation}\boldsymbol{\nabla }\boldsymbol{u}=\dot{\epsilon }\, \textit{diag}\left(-\frac{1}{2},-\frac{1}{2}\, ,1\right), \quad \dot{\boldsymbol{\gamma} }=\boldsymbol{\nabla }\boldsymbol{u}+\left(\boldsymbol{\nabla }\boldsymbol{u}\right)^{T}=\dot{\epsilon }\, \textit{diag}\left(-1,-1\, ,2\right).\end{equation}

The problem is reduced to solving the two components of the momentum balance in the $r\hbox{-}$ and $z\hbox{-}$ directions, noting that there is no spatial dependence under this idealised flow condition

(4.34) \begin{align} \frac{{\textrm{d}}\tau _{\textit{zz}}}{{\textrm{d}}t}&=2\dot{\epsilon }\tau _{\textit{zz}}+\frac{2}{Ec}\dot{\epsilon }-\frac{\max \!\left(0,\, \left| \boldsymbol{\tau} _{\boldsymbol{d}}\right| -Y_{s}\right)^{\frac{1}{n}}}{Ec}\frac{\tau _{\textit{zz}}}{\left| \boldsymbol{\tau} _{\boldsymbol{d}}\right| }, \end{align}

(4.35) \begin{align} \frac{{\textrm{d}}\binom{\tau _{\textit{rr}}}{\tau _{\theta \theta }}}{{\textrm{d}}t} &=-\dot{\epsilon }\binom{\tau _{\textit{rr}}}{\tau _{\theta \theta }}-\frac{1}{Ec}\dot{\epsilon }-\frac{\max \! \left(0,\left| \boldsymbol{\tau} _{\boldsymbol{d}}\right| -Y_{s}\right)^{\frac{1}{n}}}{Ec}\frac{\binom{\tau _{\textit{rr}}}{\tau _{\theta \theta }}}{\left| \boldsymbol{\tau} _{\boldsymbol{d}}\right| }. \end{align}

The system of the three ordinary differential equations is solved numerically. The time evolution of $\tau _{\textit{zz}}$ is shown on figure 3, for $n=0.45$ . The important finding from this simple analysis is that $\tau _{\textit{zz}}$ scales like $(t_{0}-t)^{n/n-1}$ , irrespective of other parameters, like $Ec,\, \text{ or }Y_{S}$ in extreme nonlinearity, in other words it depends only on $n$ .

Figure 3. Evolution of the axial normal stress, $\tau _{\textit{zz}}$ , vs. time to pinch in uniaxial extension flow with $\dot{\epsilon }={2}/{t_{0}-t}$ . The other EVP parameters are $n=0.45,\, Y_{S}=1$ , resulting in ${n}/{n-1}=-0.8181$ .

The same figure demonstrates that elasticity through $Ec$ affects only the prefactor of the asymptotic form, or, in other words, the intersection points of the asymptotic lines with the ordinate. Reducing $Ec$ will require longer time to reach the region of high nonlinearity, where $\tau _{\textit{zz}}$ scales as we calculated, but at the same point in time the axial stress is smaller.

Moreover, multiple tests were conducted by varying the exponent $n$ with the same extension rate $\dot{\epsilon }=({2}/{T})$ , and with $Ec=0.1$ and $Y_{s}=1$ to determine if the exponent in the scale $(t_{0}-t)^{{\alpha _{3}}}$ will agree with our prediction of $({n}/{n-1})$ . Figure 4, clearly shows that this is the case for $n\in [0.3,\, 0.8].$

Figure 4. The scaling exponent of axial normal stress, $\tau _{\textit{zz}}$ , vs. the SHB strain-rate-thinning exponent (n). The other EVP parameters are $Ec=0.1,\, Y_{S}=1$ .

5.1. Accuracy of the numerical method

To validate our numerical method, we compare our simulation results with those reported by Suryo & Basaran (Reference Suryo and Basaran2006) for a power-law material. In such a material, the viscosity is expressed as

(5.1) \begin{equation}\tilde{\mu }\left(\widetilde{\dot{\gamma }}\right)=\tilde{k}\left| \widetilde{\dot{\gamma }}\right| ^{n-1},\end{equation}

where $\tilde{\mu }$ is the apparent power-law viscosity, $\tilde{k}$ is the consistency index and $n$ is the power-law exponent. For flow without inertia, the aforementioned study predicts different flow regimes, depending on the value of $n$ . Among them the most challenging one to capture is the so-called non-slender viscous power-law regime arising for $n\lt 0.54$ , because the filament loses its slenderness and as $n$ decreases, the strongly shear-thinning behaviour of the material leads to very fast pinch-off of the thread under capillary action. In this regime, the asymptotic scales are

(5.2) \begin{align} h\, &\sim \, T^{n}, \end{align}

(5.3) \begin{align} Z\, &\sim \, T^{n}, \end{align}

(5.4) \begin{align} u_{z}\, &\sim \, T^{n-1}. \end{align}

We select $n=0.4,\, Oh^{-2}=0,\, L=2$ to validate our numerical framework. The initial interface height is set as $h(z,0)=1-0.8\cos ({\pi z}/{2})$ . This perturbed interface suffices to initiate instability, even with a stationary upper disk ( $u_{\textit{plate}}=0$ ), because of the absence of plasticity. Figure 5(a) shows the interface profile just before the breakup obtained with our numerical method and figure 5(b) illustrates a close-up view of the interface near the pinch point. Furthermore, figure 5(c) depicts the scaling of the minimum radius, axial velocity and axial position for the power-law thread. All scales in this figure are in perfect agreement with those reported in Suryo & Basaran (Reference Suryo and Basaran2006). This demonstrates the validity and accuracy of our calculations.

5.2. General time Evolution of an EVP filament

We use the 0.2 % Carbopol solution produced and rheologically characterised by Lopez et al. (Reference Lopez, Naccache, De and Mendes2017), as the material for our base case study. We used a nonlinear regression to determine the material parameters from their flow curve and frequency sweep test at 1 % strain amplitude, and we achieved close agreement with the rheological measurements. For a more detailed description, the reader is referred to Moschopoulos et al. (Reference Moschopoulos, Spyridakis, Varchanis, Dimakopoulos and Tsamopoulos2021, Reference Moschopoulos, Kouni, Psaraki, Dimakopoulos and Tsamopoulos2023). The resulting parameter values for the base case study, are given in table 2, except for setting $Oh^{-2}=0$ , to simulate flow physics without inertia.

Table 2. Parameters for the base case study.

Figure 6 illustrates the time evolution of the EVP filament. Once the disks begin to move, the material yields in most of its length, as seen for example at $t=1.32$ . Near the disks though, the material only deforms elastically as a Kelvin–Voigt solid, because its initial velocity and stresses were set to zero everywhere and the no-slip and no-penetration conditions adhere the material to the disks, preventing the development of normal stresses as well. To identify the yielded regions, we employ the von Mises criterion, as detailed in § 2.1. Yielding occurs first at the plane of symmetry, where the stress magnitude is maximised and with the assistance of capillarity they form a neck there. At all times, the filament deformation retains the plane of symmetry at $z=0$ , because of the absence of gravity and inertia. As stretching continues, the neck gets thinner there making the velocity and stress gradients more pronounced. Simultaneously, rate thinning decreases the viscous resistance there, making the deformation larger and more local, $\, t=2.75$ . Consequently, deformation has to decrease away from the mid-plane, reducing the velocity and stress gradients. Where the stress magnitude drops below the yield condition again, the material will solidify, permitting deformation in a decreasing length of the bridge. This sequence of events can be seen in the first three panels of figure 6.

Figure 6. Time evolution of filament under constant stretching with $U$ = 0.39. The red and blue areas indicate yielded and unyielded regions, respectively. The other parameters are $Oh^{-2}=0,\, n=0.45,\, Y_{s}=1.25,\, Ec=0.14,\, L_{0}=2$ .

When $h_{\textit{min}}$ reaches approximately $O(10^{-2})$ , the filament begins to deform even more rapidly, now driven primarily by capillarity. The faster escape of the material from the plane of symmetry through the increasingly smaller cross-section of the neck produces a jet-like, larger axial velocity away from the pinching region. With its increased axial velocity, this material now penetrates further inside the unyielded material away from the necking region, and toward the top moving solid disk, compressing it more effectively. In this area, the local increase in the radius is attributed to the continuous injection of material. This complex kinematics leads to material yielding away from the neck with a small unyielded area in between at time $t=2.8978$ . Before turning into the dynamics, we note the following: Because of axial symmetry, the radial velocity is zero at the axis of symmetry, but away from it and near the pinch point, at $z=0$ , it is negative following the decrease in radius, whereas in the corner region (close to highest curvature of the interface), it becomes positive. Consequently, close to the pinch point, $\tau _{\textit{rr}}$ is negative but it turns positive in the vicinity of the corner region, explaining the sign change of $\tau _{\textit{rr}}$ in this area; see figure 7 below. The same reasoning applies to the change of the sign of $\tau _{\textit{zz}}$ and $\tau _{rz}$ . The sign change of all stresses indicates that, close to the corner region, the components of the deviatoric stress tensor will approach zero. Hence, the transition region contains unyielded material that separates the yielded neck from the yielded island further up. As the filament approaches breakup, velocity and stress magnitude increase even faster and the unyielded area in the transition region shrinks in size and eventually disappears, as shown in figure 6, at $t=2.8984$ . It is worth mentioning that this is observable in EVP materials, because they can deform even in unyielded regions.

Figure 7. Contours of velocity and stress components of the filament under constant stretching with $U$ = 0.39 just before pinch-off, at $t=2.8984$ . Other parameters are $Oh^{-2}=0, n=0.45, Y_{s}=1.25, Ec=0.14, L_{0}=2$ .

Figures 7(a) and 8(a) present contours of velocity and extra stress tensor components, along with the second invariant of the extra stress tensor, just before pinch-off at $t=2.8984$ . These plots clearly illustrate that extremely high velocity and stress magnitudes arise in the necking region. Moving away from this region, all values rapidly decrease until they match the boundary conditions imposed on the solid plate. Figures 7(b) and 7(c) depict contours of velocity components focusing at the corner and necking regions of the interface, respectively. In the necking region, the axial velocity $u_{z}$ dominates, becoming three orders of magnitude larger than the radial velocity $u_{r}$ . This dominance decreases when approaching the corner region, with its magnitude being only one order of magnitude higher than that of $u_{r}$ .

Figure 8. Contours of pressure and stress components of the filament under constant stretching with $U$ = 0.39 just before pinch-off, at $t=2.8984$ . Other parameters are $Oh^{-2}=0, n=0.45, Y_{s}=1.25, Ec=0.14, L_{0}=2$ .

As demonstrated in figures 7(c) and 7(d), the magnitudes of $\tau _{\textit{rr}}$ and $\tau _{\theta \theta }$ are nearly identical in both magnified regions, with their extrema occurring at the necking region. This is consistent with the fact that $({{\textrm{d}}u_{r}}/{{\textrm{d}}r})$ and $({u_{r}}/{r})$ also attain their maximum value in this region. The pressure magnitude in the necking region is proportional to $({1}/{h})$ , implying that the high pressure within this region decreases progressively as the distance from the symmetry plane increases. This trend is clearly illustrated in figures 8(b) and 8(c), which show pressure contours in the magnified corner and necking regions, respectively. The shear stress component $\tau _{rz}$ reaches its extremum in the corner region, approximately at the point of highest interface curvature. Figure 8(b) further reveals that $\tau _{rz}$ exhibits a magnitude comparable or even slightly higher than $\tau _{\textit{rr}}$ and $\tau _{\theta \theta }$ in this region. The extension-dominated nature of filament breakup is evident in the contours of $\tau _{\textit{zz}}$ in necking and corner regions, as shown in figure 8(c) and 8(d). Similar to the axial velocity, $\tau _{\textit{zz}}$ is three orders of magnitude higher than all other stress components in the necking region but rapidly decreases until it becomes comparable to $\tau _{rz}$ at the corner. The magnitude of the second invariant of the extra stress tensor in the thread is approximately $({1}/{\sqrt{3}})| \tau _{\textit{zz}}|$ , see figure 8(c); However, this relation no longer holds beyond the neck region, as depicted in figure 8(d). The variation in the radial direction inside the neck for all variables is negligible. One may observe a small radial variation of $u_{r}, \tau _{rz }\text{ and}\, \tau _{\textit{rr}}$ , but only because their magnitude is much smaller than that of the other variables and the range of colours in each panel depends on the range of variation of the corresponding variable. In the supplementary material SM2, we present the time evolution of the axial stress component $\tau _{\textit{zz}}$ and the axial velocity $u_{z}$ , where one can easily observe that they increase from $O(1)$ to $O(10^{5})$ without any oscillations in their contours, during the same simulation.

All these justify the assumptions made to derive the slender filament equations in the necking region. However, we will not employ the slender filament equations for simulating this flow. As stated by Moschopoulos et al. (Reference Moschopoulos, Syrakos, Dimakopoulos and Tsamopoulos2020), in stretching of viscoplastic filament, the deviation between 2-D and slender filament simulations increases as the viscoplastic character of the flow becomes more pronounced. We believe this applies for EVP material as well and verifying the accuracy of the solution obtained by the slender equations requires extensive analysis, particularly for different combinations of $Y_{s}$ and $Ec$ . Moreover, aside from the necking region, the filament is not slender throughout, particularly in the corner region, necessitating a full 2-D analysis for accurate results. Additionally, the slender filament approximation, neglects $\tau _{rz}$ , which is inappropriate given that $\tau _{rz}$ is of the same order of magnitude or larger than $\tau _{\textit{rr}}$ and $\tau _{\theta \theta }$ in this region; see also Eggers, Herrada & Snoeijer (Reference Eggers, Herrada and Snoeijer2020).

5.3. Determination of pinch-off time, ${t}_{{0}}$

Accurate approximation of the pinch-off time, $t_{0}$ , is most important, because even a very small error in determining its value will lead to incorrect calculation of the asymptotes. Moreover, in our arbitrary Lagrangian Eulerian formulation, breakup cannot occur without external intervention, which must take place as close as numerically possible and physically allowed to the pinch-off time. Computational experiments as well as other studies suggest that, in order to accurately determine it, simulations have to be conducted until the minimum radius becomes approximately five orders of magnitude smaller than the initial radius, or, in the present dimensionless form, when $h_{\textit{min}}=10^{-5}$ . Then the very last time instant of the simulation can be considered as pinch-off time and the scalings are correctly determined. On the other hand, when the minimum radius reaches approximately $O(10^{-3})$ , the EVP stresses grow very fast because of the viscoelastic nature of the EVP material, making the resulting system of equations extremely stiff to solve. Although reaching values as small as $h_{\textit{min}}=10^{-5}$ may not be feasible, we have employed a least-squares, nonlinear regression method to extrapolate the pinch-off time.

The extrapolation to the breakup time is based on the asymptotic form obtained in § 4.1 for $h_{\textit{min}}=\phi _{0}(t_{0}-t)^{{\alpha _{1}}}$ . Treating $\phi _{0}, t_{0}$ and $\alpha _{1}$ as fitting parameters, we apply this power-law function to the minimum radius data over the last 100 computed time steps. Figure 9 depicts the regression applied to the minimum radius as a function of time for the base case study. The computed $t_{0}$ for this case is found to be $2.8984063$ , while $\varPhi _{o }\approx 10.7$ and $\alpha _{1}=0.8$

Figure 9. Nonlinear regression employed to extract the pinch-off time, $t_{0}$ , from the numerical results of the base case study, for $Oh^{-2}=0, n=0.45, Y_{s}=1.25, Ec=0.14, L_{0}=2$ .

The calculated $t_{0}$ is then used to generate plots of $h_{\textit{min}}, z, u_{z}, \tau _{\textit{zz}}$ versus $T$ . According to § 4.1, the axial stress $\tau _{\textit{zz}}$ should scale as $(t_{0}-t)^{({n}/{n-1})}$ . Owing to conservation of momentum and the balance between capillarity and axial stresses, the minimum radius should decrease asymptotically as $({1}/{\tau _{\textit{zz}}})$ , in other words as $h_{min }\sim (t_{0}-t)^{({n}/{1-n})}.$ Furthermore, to satisfy the slender jet continuity equation, the axial position and velocity must scale like $(t_{p}-t)^{\delta }$ and $(t_{p}-t)^{1-\delta }$ , respectively. Figure 10 illustrates in log–log scales $h_{\textit{min}}, z, u_{z}, {\textrm{and}}\, \tau _{\textit{zz}}$ versus $t_{0}-t$ for $n=0.45, Ec=0.14, Y_{s}=1.25, U=0.39$ . The linear slope in figure 10 verifies that, indeed, a power-law relation exists for all variables. The calculated slopes are $h_{\textit{min}}, Z, u_{z}\, {\textrm{and}}\, \tau _{\textit{zz}}$ are $0.80\, (\approx ({n}/{1-n})),-0.08\, (=\delta ), -1.08\, (=\delta -1)$ and $-0.81\, (\approx {n}/{n-1})$ , respectively. These slopes represent the exponents of power-law scaling functions for the base material and they are in perfect agreement with the asymptotes we derived with the elasto-capillary scalings, as discussed earlier.

Figure 10. Decay of the minimum radius, growth of axial velocity, axial position and axial stresses as pinching time is approached for $Oh^{-2}=0, n=0.45, Ec=0.14, Y_{s}=1.25, U=0.39$ ( $({n}/{1-n})=0.818$ . The slopes determine the corresponding exponents of the power-law scaling functions.

5.4. Self-similar forms of the two-dimensional solutions

Given a distinct set of parameters $(n,Ec,Y_{s})$ , first we calculate the scaling exponents for minimum radius, axial position, axial velocity and stress, as shown in figure 10, then we rescale these variables accordingly at different times to verify their collapse onto a single curve, completing in this way the demonstration of the existence of self-similarity. We have observed that this holds for a range of power-law exponents, but, for conciseness, next we present only a single set of parameters corresponding to the base case study. The inset to figure 11 shows interface profiles for five different values of minimum radius $h_{\textit{min}}$ . Figure 11 depicts that these transient interface profiles obtained by 2-D simulations collapsing onto a single self-similar interface profile by rescaling the interface height $\overline{h}=({h}/{(t_{0}-t)^{({n}/{1-n})}})$ and the axial position $\overline{Z}=({z-z_{0}}/{(t_{0}-t)^{\delta }})$ , where $\delta \approx -0.08$ . Please note that the rescaled variables are denoted by an overbar ( $\overline{})$ over the corresponding variable symbol.

Figure 11. (Inset) Transient interface shape obtained by 2-D simulations at five different $h_{\textit{min}}$ values. (Main figure) Rescaled interface using the numerically determined scaling laws for the same five $h_{\textit{min}}$ values. The main figure illustrates convergence to a self-similar interface profile. The very large difference in the magnitudes of $h$ before and after scaling is noteworthy.

The inset to figure 12 shows the variation of axial velocity $u_{z}$ at the interface versus the axial coordinate $z$ for five different values of minimum radius $h_{\textit{min}}$ . Figure 12 depicts that by rescaling the axial velocity as $\overline{u}_{z}= ({u_{z}}/{(t_{0}-t)^{{\alpha _{2}}}})$ , with $\alpha _{2}=-1.08$ , the transient velocity profiles collapse onto a single self-similar profile as the minimum radius $h_{\textit{min}}\rightarrow 0$ .

Figure 12. (Inset) Transient velocity variation along axial coordinate (z) at the interface obtained by 2-D simulations at five different $h_{\textit{min}}$ values. (Main figure) Rescaled velocity plotted against rescaled axial coordinate using the numerically determined scaling laws for the same five $h_{\textit{min}}$ values. The main figure demonstrates convergence to a self-similar velocity profile. The very large difference in the magnitudes of the velocity before and after scaling is noteworthy.

Finally, the inset to figure 13 shows the axial stress $\tau _{\textit{zz}}$ versus the axial coordinate extracted from nodal values of 2-D calculations at the interface for five different values of $h_{\textit{min}}$ . The rescaled axial stress $\overline{\tau }_{\textit{zz}}=({\tau _{\textit{zz}}}/{(t_{0}-t)^{({n}/{n-1})}})$ vs. rescaled axial coordinate does converge to a self-similar profile, as shown on figure 13.

Figure 13. (Inset) Transient axial stress variation along axial coordinate (z) at the interface obtained by 2-D simulations at five different $h_{\textit{min}}$ values. (Main figure) Rescaled axial stress plotted against rescaled axial coordinate using the numerically determined scaling laws for the same five $h_{\textit{min}}$ values. The main figure demonstrates convergence to a self-similar axial stress profile. The very large difference in the magnitudes of stress before and after scaling is noteworthy.

The plots in this section show good collapse near $\overline{Z}=0$ ; however, the quality of the collapse diminishes as $\overline{Z}$ increases. This can be attributed to the negative value of the exponent $\delta$ , which leads to smaller values of $\overline{Z}$ compared with cases where $\delta \gt 0$ . As illustrated in the scaling plot of figure 10, the axial position increases as breakup is approached. Consequently, in self-similarity plots, the span of $\overline{Z}$ is reduced relative to other constitutive equations, such as Newtonian or power-law fluids, for which $\delta \gt 0$ .

Moreover, we have extracted numerically scaling laws for the variables in the corner region. However, they do not follow a universal scaling law, we were unable to derive analytically these scalings nor could we obtain self-similar solutions. It is possible that this is caused by the fact that the corner region involves in the thread side yielded material until the very late stages of breakup, much later than when asymptotic behaviour of the neck starts, but unyielded material in the other side toward the disc. Consequently, derivation of the asymptotics is not straightforward, because the material exhibits dual behaviour in this region. Additionally, since the solvent viscosity is set to zero in this study, the dominant balance and self-similar form of the governing equations become more complex than those reported by Fontelos & Li (Reference Fontelos and Li2004) for the Giesekus fluid.

5.5. Parametric analysis

We will separate this section based on the value of the power-law index. When $0.43\lt n\lt 1$ , surface tension is the main driving force causing breakup. This range of values is examined in § 5.5.1. Although for $n\gt 0.5$ , the SHB model predicts extension rate hardening, which is not expected for Carbopol (Kordalis et al. Reference Kordalis, Varchanis, Ioannou, Dimakopoulos and Tsamopoulos2021), we will examine values of $n$ up to one for a more complete study. Values of $0.33\lt n\lt 0.43$ will be examined separately in § 5.5.2, because then capillarity plays no role in pinching-off, while $n=1$ will be examined in § 5.5.3, also separately, because, as we will demonstrate, pinching-off does not take place in the Saramito–Bingham model (Saramito Reference Saramito2007). We will examine the effect of $n$ and $Ec$ , where appropriate. Here, $Y_{s}$ has no effect, because the thread during pinching-off is yielded throughout.

5.5.1. Elasto-plasto-capillary breakup regime, ${0.43}\lt {n}\lt {1}$

Although the YS does not enter the dominant balance determining the asymptotic behaviour of the thread in this regime, the title includes the term ‘plasto’, because part of the included elastic term is what remains from the viscoplastic contribution in the constitutive law, see (3.11).

5.5.1.1. Effect of strain-rate-thinning exponent (n)

We proceed by examining the effect of the strain-rate-thinning exponent, $n$ , on the pinch-off dynamics. Figure 14(a) depicts interface profiles for four different values of $n$ , recorded when $h_{\textit{min}}=1.5\times 10^{-3}$ .

Figure 14. Profiles of the entire filament (a) and closeup views in the corner (b) and necking (c) regions for four different values of $n$ . In (b) we did not include the entire ‘local’ interface for $n=0.8$ , because it would reduce the clarity of the rest. The remaining parameter values are $Y_{s}=0.1, Ec=0.14, U=0.1$ .

As $n$ increases, strain-rate thinning decreases, resulting in increasing elastic stresses, including its radial normal component, which opposes capillarity. This leads to a longer thread around the middle of the filament for the same $h_{\textit{min}}$ , as depicted more clearly in figure 14(b). Given that the initial dimensions of the filament for all cases are the same, a longer neck causes a slight increase in the total height of the filament (see figure 14 a), when it has reached the same $h_{\textit{min}}$ under the same constant stretching velocity. Additionally, the increased opposing elastic force slows down the breakup rate, causing the EVP thread to break later in time as $n$ increases. Figure 14(c) shows that the interface profile close to pinch off is almost a perfect cylinder for all values of $n$ .

Figure 15. Effect of the exponent of strain-rate thinning on the scaling with respect to $(t_{o}-t)$ of minimum radius (a), axial stress (b), axial velocity (c) and axial length of neck (d). The rest of the parameter values are $Y_{s}=0.1, Ec=0.14, U=0.1$ .

The dependence on time to pinch of the minimum radius, the axial normal stress and axial velocity, as well as the length of the neck up to which its radius increases by only 1 % from its minimum value, are presented in figure 15. Figures 15(a) and 15(b) depict the $h_{min }$ decrease and the $\tau _{\textit{zz}}$ increase as a function of time to pinch $t_{o}-t$ . Focusing on their values when they reach their asymptotic behaviour, one readily observes that they follow very closely the power-law dependence on $n$ , namely $\pm n/(1-n)$ for all four values of $n$ . The calculated slopes, which represent the exponents of the power-law scaling functions, are almost in perfect agreement with the proposed asymptotes. An error of $\pm \, 2\,\%$ from the earlier predicted asymptotes may be attributed to inaccuracies of the numerical method and a slight radial dependence of the variables in the numerical solution of the 2-D equations. As these exponents require, when $n$ increases, the slopes of the asymptotes of both variables increase. Figure 15(a) also shows that increasing $n$ brings $h_{min }$ faster to a certain small value, although the time to breakup is delayed, since this occurs at larger values of $(t_{o}-t)$ . Similarly, figure 15(c), demonstrates that the slope of the asymptotic increase of $v_{z}$ increases with $n$ , although for $n\lt 0.5$ , the exponent does not have an analytic relation with it; see (3.27) and (3.28). Finally, figure 15(d) indicates that the thread of the filament becomes longer at a faster rate, when $n$ increases. This length extends from the plane of symmetry to the axial position where $h=1.01h_{min }$ , and it is denoted as $Z$ .

Altogether, figure 15 demonstrates that, for smaller $n$ values, the asymptotic behaviour is reached when calculations are carried out to times much closer to breakup. This leads to a thinner neck, driven by higher velocity produced by higher axial stress. Another important observation from these four panels is the sharp increase in the extreme values reached by $h_{\textit{min}}, v_{z}\, {\textrm{and}}\, \tau _{\textit{zz}}$ as $n$ decreases. These values are $\sim 2\times 10^{-4}, 2\times 10^{4}\, \text{and }5\times 10^{4}$ , respectively, for $n=0.45$ , while their initial values are 1, 0 and 0, respectively. This makes plain how challenging these calculations are and the importance of the new ideas adopted in the particular numerical algorithm.

5.5.1.2. Effect of elastic modulus

We move on by examining the effect of the elasto-capillary number, $Ec$ , on the pinch-off dynamics. Figure 16 illustrates the different shapes of the filament when $h_{\textit{min}}=4.5\times 10^{-4}$ .

Figure 16. Profiles of the entire filament (a) and closeup views in the corn (b) and necking (c) regions for three different values of $\, Ec$ . The rest of the parameter values are $Y_{s}=0.1, n=0.45, U=0.1$ .

Focusing on figure 16(a), we observe that increasing $Ec$ generates both a shorter filament and, hence, one that breaks faster, but a longer thread around the plane of symmetry. This seemingly contradictory trend is explained as follows: increasing $Ec$ affects both the unyielded and yielded regimes of the filament in different ways. To start with, when $Ec$ increases, the filament deforms more as an elastic solid prior to yielding, allowing for a more localised viscoelastic deformation in the region around the plane of symmetry. This leads to earlier formation of the neck and earlier breakup of the filament, although the opposing elastic forces also increase. All these lead to a shorter total filament height at pinch-off time, figure 16(a). This increased elastic yielding has been reported also by Moschopoulos et al. (Reference Moschopoulos, Kouni, Psaraki, Dimakopoulos and Tsamopoulos2023), where the bulk dynamics has been studied. The same increase in $Ec$ increases $\tau _{\textit{zz}}$ , which accelerates the evacuation of the neck and the decrease of its radius. At the same time, the larger normal stresses generate a longer neck, because they oppose neck thinning by capillarity, as shown on closeup view at the corner region, figure 16(b). Figure 16(c) shows that the neck profile close to pinch off forms a nearly perfect cylinder for all values of $Ec$ .

The trends in $h_{\textit{min}}$ decay and $\tau _{\textit{zz}}, v_{z},\text{ and}\, Z$ growth are presented in figure 17. The slopes of all four variables with respect to the time to pinch $(t_{o}-t)$ in the asymptotic regime depend only on $n$ , but are independent from $Ec$ . The elastocapillary number affects only the intersection with the ordinate, via affecting the coefficient $\phi _{0}$ in (3.20). Even this effect is weak for all variables, except for $Z$ . Indeed, figure 17(d) indicates that the axial extent over which this scaling holds is longer, or, in other words, the neck is longer for higher $Ec$ , although its growth rate remains the same. The calculated slopes, which represent exponents of the power-law scaling functions, are almost in perfect agreement with the proposed asymptotes.

Figure 17. Effect of $Ec$ on the scaling with respect to $(t_{o}-t)$ of minimum radius (a), axial stress (b), axial velocity(c) and axial length of neck (d). The rest of the parameter values are $Y_{s}=0.1, n=0.45, U=0.1$ .

5.5.2. Elasto-plastic breakup regime, ${0.33}\lt {n}\lt {0.43}$

In this section we will examine the pinching dynamics of materials with smaller power-law exponent $0.33\lt n\lt 0.43$ and with as low $n$ values (always positive) as our algorithm allows with convergence in a reasonable time. As in the previous regime, the YS does not enter the dominant balance in this regime either, but its title again includes the term ‘plastic’, because part of the appearing elastic term is what remains from the viscoplastic term in the constitutive law. A small decrease in $n$ below 0.45 leads to a significant deviation from the asymptotic result for the minimum radius. For the values of $0.33\lt n\lt 0.43$ , the breakup mechanism is no longer governed by capillarity. Instead, it is dictated purely by elastic stresses, which drive the breakup after the neck is formed. This type of pinching-off has been identified by Renardy (Reference Renardy2002a ) for the generalised form of PTT materials that he introduced. To further investigate this regime, we carried out filament stretching simulations for three different values of $n$ where $0.33\lt n\lt 0.43.$ In the elasto-plastic breakup regime, the minimum radius follows a different asymptotic scaling $h_{\textit{min}}\sim \, (t_{0}-t)^{({n}/{4-8n})}$ . However, since surface tension does not appear in the constitutive equation, the scaling for $\tau _{\textit{zz}}$ remains unchanged, following $\tau _{\textit{zz}}\sim (t_{0}-t)^{({n}/{n-1})}$ .

Figure 18. Effect of the exponent of the strain-rate thinning (n) in elastic breakup regime on the scaling with respect to $t_{0}-t$ of minimum radius (a), axial stress (b), axial velocity (c) and axial length of neck (d). The rest of the parameter values are $Y_{s}=0.1, Ec=0.14, U=0.1$ .

This is shown in figures 18(a) and 18(b), respectively. Additionally, our calculations indicate that, by decreasing the value of $n$ , the axial velocity and axial position scaling exponents, $\alpha _{2}$ and $\delta$ approach the values of −1 and 0, respectively, as illustrated in figures 18(c) and 18(d), respectively. The exponents in this regime do not depend on $Ec$ , similar to the capillary-driven regime, and the corresponding results are omitted for conciseness. It is noteworthy that reaching the asymptotic range now requires getting even closer to the breakup time $(t_{0}-t)\lt 10^{-5}$ and fully capturing it requires approaching small $T$ values more than ever before, down to $(t_{0}-t)\sim 10^{-8}$ . All these necessitate extremely small time steps for accurate predictions. Furthermore, the extreme values of $h_{\textit{min}}, v_{z}\, {\textrm{and}}\, \tau _{\textit{zz}}$ now become $\sim 3\times 10^{-4}, 5\times 10^{5}\, \text{and }2\times 10^{5}$ , respectively, with the last two approximately one order of magnitude larger than their values for $n=0.45$ . All these observations reveal that the calculations become increasingly difficult as $n$ decreases. The deviation in the scaling exponent $\alpha _{1}$ for $n=0.4$ is most pronounced among all cases studied. Specifically, our numerical results yield $\alpha _{1} = 0.45$ , whereas the asymptotic scale predicted from slender filament equations is 0.5. Even with a more refined mesh in the necking region, including more than 10 000 axial elements results in the same value of $\alpha _{1}=0.45$ . We believe this deviation arises by the factors discussed earlier.

5.5.3. Collapse of filament shapes and radial velocity asymptotes for all cases with, ${n}\lt {1}$

Having analysed the scalings where the SHB exponent, $n$ , is less than unity, we now present self-similar interface profiles obtained by 2-D simulations for various values of $n$ , shown in figure 19. We have shown that $z$ scales as $(t_{0}-t)^{\delta }$ and that when $n$ increases, $\delta $ increases in absolute value, remaining negative. Consequently, the rescaled axial coordinate $\overline{Z}$ spans a narrower region as $n$ increases. Similarly, the span of $\overline{h}$ decreases. Both are clearly seen in figure 19. The local slenderness of the necking region decreases progressively with decreasing $n$ . This trend can be understood by examining the ratio of radial to axial length scale, given by $({h_{\textit{min}}}/{z})\sim (t_{0}-t)^{{\alpha _{1}}-\delta }$ . The variation of the exponent $\alpha _{1}-\delta$ is provided in table 3, revealing a sharp decrease from ∼ 7 to 0.45 as $n$ decreases from $0.8$ to $0.4$ . Consequently, the breakup dynamics in the necking region exhibits an abrupt transition from a fully acceptable slender filament one to one that increasingly approaches a non-slender breakup at $n=0.4$ .

Table 3. Values of exponents $\alpha _{1}, \delta$ and $\alpha -\delta$ for five different SHB exponents. The rest of parameters are $Y_{s}=0.1, Ec=0.14, U=0.1$ .

Figure 19. Self-similar interface profiles of EVP fluids with five different strain-rate-thinning exponents. The rest of the parameters are $Y_{s}=0.1, Ec=0.14, U=0.1$ .

Figure 20 shows the evolution of the radial velocity at the interface on the symmetry plane as a function of remaining time to pinch, $t_{0}-t$ . Results from 2-D simulations reveal that, for all values of the SHB exponent $(n)$ less than unity, the scale of $u_{r}$ is $\alpha _{1}-1$ . So, the radial velocity exhibits an asymptotic scale $u_{r}\, \sim ({h_{\textit{min}}}/{(t_{0}-t}))$ , being valid for all cases where $n\lt 1$ . Interestingly, as shown in figure 20, for SHB exponents $n\gt 0.5$ , the absolute radial velocity $| u_{r}|$ at the interface on the symmetry plane decreases considerably as $t_{0}-t\rightarrow 0$ , indicating an increase in material’s resistance to breakup, which is a characteristic of extension-rate-thickening behaviour. Conversely, for $n\lt 0.5$ , $| u_{r}|$ increases toward breakup, reflecting a significant reduction in elastic resistance, corresponding to extension-rate-thinning behaviour.

Figure 20. Variation of the radial velocity magnitude at the intersection between the symmetry plane and the interface with $t_{0}-t$ . The rest of the parameters are $Y_{s}=0.1, Ec=0.14, U=0.1$ .

5.5.4. Discrimination between finite-length and point breakup

In this section we elaborate further on the discussion in § 4.1 concerning the finite-length breakup. We investigate this phenomenon by examining the filament thickness $h$ at two axial locations; one at the plane of symmetry, where $h$ attains its minimum value, and another one slightly above it at $z^{+}=0.01$ . Initially, we attempted to demonstrate finite-length breakup by plotting $h_{z+}$ and $h_{\textit{min}}$ as functions of time $t$ or time to breakup $t_{0}-t$ . However, the near perfect overlap of these curves hindered any definitive conclusion. Next, we employed an alternative metric by plotting the normalised difference of filament thickness probed at these two axial locations versus time to breakup

(5.5) \begin{equation}\chi _{h}=\frac{h_{{z^{+}}}-h_{\textit{min}}}{h_{\textit{min}}}.\end{equation}

This quantity may also be experimentally accessible. If a line is drawn connecting these two points on the interface, one at $z=0$ and the other at $z=z^{+}=0.01$ , the difference $h_{{z^{+}}}-h_{\textit{min}}$ is proportional to the absolute slope of that line, since ${\unicode[Arial]{x0394}} z$ is fixed. A decrease in $h_{{z^{+}}}-h_{\textit{min}}$ indicates that the filament shape is becoming more cylindrical. Meanwhile, the decay rate of $h_{\textit{min}}$ characterises the rate at which breakup occurs at the symmetry plane. Therefore, an increase in $\chi _{h}$ as $t_{0}-t\rightarrow 0$ denotes that the rate of filament thinning leading to point breakup exceeds the rate at which a cylindrical ‘inner thread’, if present, is formed. In this context, a negative slope of $\chi _{h}$ plotted against $t_{0}-t\rightarrow 0$ in log–log plot implies that breakup will asymptotically precede the formation of a cylindrical neck, if such neck is formed at all. Conversely, a positive or zero slope suggests finite-length breakup, wherein the filament retains a finite axial extent close to breakup. Asymptotically, this implies that the cylindrical neck forms either before breakup, if $\chi _{h}\rightarrow$ 0, or simultaneously with it, if $\chi _{h}$ approaches a constant. As shown in figure 21, $\chi _{h}$ increases near breakup for both a power-law material with $n=0.4$ analysed in § 5.1 and a Newtonian fluid. Conversely, in an EVP filament with $n$ varying from 0.4 to 0.6, analysed in §§ 5.5.1 and 5.5.2, $\chi _{h}$ exhibits a decreasing or constant trend, with a plateau at O( $10^{-2}-10^{-4})$ as breakup is approached. These observations suggest that the EVP filament retains an almost cylindrical shape very near pinch-off, consistent with finite-length breakup behaviour.

Figure 21. Evolution of $\chi _{h}=({h_{{z^{+}}}-h_{\textit{min}}})/{h_{\textit{min}}}$ plotted against time to breakup $t_{0}-t$ for EVP, power-law and Newtonian fluids. Other parameters for the EVP fluid are $Y_{s}=0.1, Ec=0.14, U=0.1$ , and those for the power-law fluid are provided in § 5.1. Here, $z^{+}=0.01$ .

5.5.5. Analysis of the local curvature at the symmetry plane

In addition to analysing the scaling laws for minimum radius, axial stress or axial velocity at the interface on the symmetry plane $z=0$ , it is often important to examine the interfacial curvature and its asymptotic behaviour at this location. This becomes particularly relevant when it is important to distinguish between a SHB material with $n=0.5$ and a Newtonian fluid, because both exhibit a similar linear decay of the minimum radius. As highlighted by Du et al. (Reference Du, Ohtani, Ellwood and McKinley2024), evaluation of the curvature at the mid-plane offers a more comprehensive characterisation of the filament shape than the decay of the minimum radius alone. Following this approach, we define the curvature ratio $\chi$ at the interface on the symmetry plane as

(5.6) \begin{equation}\chi = \displaystyle\frac{\kappa _{z}}{\kappa _{r}}= \displaystyle\frac{\left| \displaystyle\frac{h_{\textit{zz}}}{\left(1+h_{z}^{2}\right)^{\frac{3}{2}}}\right| }{\left| \displaystyle\frac{1}{h\left(1+h_{z}^{2}\right)^{\frac{1}{2}}}\right| }.\end{equation}

Substituting (3.3) and (3.6) into the above expression yields a general asymptotic scaling form for $\chi$

(5.7) \begin{equation}\chi \sim \left(t_{0}-t\right)^{{2\alpha _{1}}-2\delta \, }\sim \, \left(h_{\textit{min}}\right)^{2-\frac{2\delta }{\alpha _{1}}}.\end{equation}

Figure 22 shows the evolution of the curvature ratio, $\chi$ , as a function of the of minimum radius $h_{\textit{min}}$ , neglecting inertial effects. It includes SHB fluids in the range of exponents $n\in [0.4,0.6]$ , as well a power-law fluid with $n=0.4$ and a Newtonian fluid.

Figure 22. Evolution of the curvature ratio against minimum radius for EVP, power-law and Newtonian fluids. Other parameters for EVP fluid are $Y_{s}=0.1, Ec=0.14, U=0.1$ , and those for the power-law fluid are provided in § 5.1.

According to Renardy (Reference Renardy1994), the inertialess breakup asymptotic scaling exponents of a Newtonian fluid are $\alpha _{1}=1\, \text{and }\delta =0.175$ , while for a power-law fluid with $n=0.4$ , they are $\alpha _{1}=\delta =0.4$ , as discussed in § 5.1. Consequently, the asymptotic scaling of $\chi$ is expected to follow $(h_{\textit{min}})^{1.65}$ for Newtonian fluid and remain constant for power-law fluids under the assumption of negligible inertia. Notably, figure 22 reveals that the EVP fluid with $n=0.5$ exhibits a numerically determined scaling exponent of approximately 2.35, which is significantly higher than that of the Newtonian case despite the fact that both fluids exhibit a similar linear decay in the minimum radius. This demonstrates that the evolution of the curvature ratio with respect to the minimum radius can effectively distinguish between the two cases, even though the formation of an elastic neck in the EVP fluids may not be fully observable in experiments. It is noteworthy that, having computed the curvature numerically, to evaluate curvature ratio $\chi$ , which involves the second derivative of the interface shape, we have calculated the $\phi _{2}$ term in (3.18) as well. We opted not to report it explicitly for conciseness.

5.5.6. Elasto-plasto-capillary regime with no breakup in finite time, ${n}={1}$

Increasing $n$ to unity to examine numerically the asymptotics of the Saramito–Bingham model, we revert to the asymptotes of the simpler Oldroyd-B fluids, which are known to predict that all variables and in particular the minimum radius have an exponential dependence on time to breakup. This is demonstrated in figure 23, where the minimum radius depends exponentially on time t, and the coefficient of the exponent is the same as the one for an Oldroyd-B fluid, $-1/(3\times Ec\times \ln (10))$ , while the coefficient of $\tau _{\textit{zz}}$ has the opposite sign. The other two variables also have an exponential dependence on $t$ with the same coefficient between them, as required by the asymptotics. The final values of all variables required to cover the asymptotic range are much smaller now than in cases with $n\lt 1$ , revealing that indeed the previous calculations are the most challenging ones. Once again, the YS plays no apparent role in the asymptotes.

Figure 23. Effect of $Ec$ on the scaling with respect to time $(t)$ of minimum radius (a), axial stress (b), axial velocity (c) and axial length of neck (d). The rest of the parameters are $Y_{s}=0.1, n=1., U=0.1$ .