4.1. Linear instability analysis based on the 1-D model

In this subsection, a brief linear analysis is performed based on the present 1-D model. Following the normal mode method, the quantities are decomposed into a base and a perturbation part, i.e.

(4.3)\begin{align} (S, u, \tau_{zz}, \tau_{rr}, A_{zz}, A_{rr}, \varGamma)&=\left(1, 0, 0, 0, \frac{L^2}{L^2+3},\frac{L^2}{L^2+3}, \varGamma_0\right)\nonumber\\ &\quad +(\hat{S}, \hat{u}, \hat{\tau}_{zz}, \hat{\tau}_{rr}, \hat{A}_{zz}, \hat{A}_{rr}, \hat{\varGamma})\exp(\alpha t+jkz)+{\rm c.c.}, \end{align}

where the hat denotes the initial amplitudes of the perturbation parts, $\alpha$ is the complex frequency with the real part $\alpha _r$ the temporal growth rate and the imaginary part $\alpha _i$ the speed of wave propagation, $j$ is the imaginary unit and ‘c.c.’ denotes the complex conjugate.

Linearizing the 1-D equations (3.1)–(3.7) and then substituting the decomposition (4.3) into the linearized equations, we obtain a set of linear homogeneous algebraic equations for the unknowns $(\hat {S}, \hat {u}, \hat {\tau }_{zz}, \hat {\tau }_{rr}, \hat {A}_{zz}, \hat {A}_{rr}, \hat {\varGamma })$. The determinant of the coefficient matrix of the equations has to be zero to ensure the system has non-trivial solutions, which yields the following dispersion relation:

(4.4)\begin{align} \alpha+3\beta Ohk^2&={-}k^2(1-\beta)Oh\frac{L^2+3}{L^2} (\varXi_1-\varXi_2)-\frac{k^2E\varGamma_0}{2(1-\varGamma_0) \left(\alpha+\dfrac{k^2}{Pe}\right)}\nonumber\\ &\quad+\frac{\sigma_ik^2(1-k^2)}{2\alpha}-\frac{B_{\kappa\infty} +9B_{\mu\infty}}{2}\varGamma_0Ohk^2, \end{align}

where

(4.5)$$\begin{gather} \varXi_1=\dfrac{\left| \begin{array}{cc} \dfrac{2L^2}{L^2+3} & \dfrac{2(L^2+3)}{L^4} \\ -\dfrac{L^2}{L^2+3} & De \alpha +\dfrac{L^2+3}{L^2}+\dfrac{2(L^2+3)}{L^4} \\ \end{array} \right| }{\left| \begin{array}{cc} De \alpha +\dfrac{L^2+3}{L^2}+\dfrac{L^2+3}{L^4} & \dfrac{2(L^2+3)}{L^4} \\ \dfrac{L^2+3}{L^4} & De \alpha +\dfrac{L^2+3}{L^2}+\dfrac{2(L^2+3)}{L^4} \end{array} \right|}, \end{gather}$$

(4.6)$$\begin{gather}\varXi_2=\dfrac{\left| \begin{array}{cc} De \alpha +\dfrac{L^2+3}{L^2}+\dfrac{L^2+3}{L^4} & \dfrac{2L^2}{L^2+3} \\ \dfrac{L^2+3}{L^4} & -\dfrac{L^2}{L^2+3} \\ \end{array} \right| }{\left| \begin{array}{cc} De \alpha +\dfrac{L^2+3}{L^2}+\dfrac{L^2+3}{L^4} & \dfrac{2(L^2+3)}{L^4} \\ \dfrac{L^2+3}{L^4} & De \alpha +\dfrac{L^2+3}{L^2}+\dfrac{2(L^2+3)}{L^4} \end{array} \right|}, \end{gather}$$

and $\sigma _i(=1+E\ln (1-\varGamma _0))$ is the initial surface tension corresponding to the initial surfactant concentration $\varGamma _0$. Note that the four terms on the right-hand side of (4.4) are related to the viscoelasticity of the polymer solution, the Marangoni stress, the surface tension and the surface shear and dilatational viscosities, respectively.

When the finite extensibility $L$ approaches infinity, the dispersion relation (4.4) reduces to that for the linear Jeffreys model and any nonlinear viscoelastic model (e.g. Oldroyd-B and FENE-P) that becomes the Jeffreys model after linearization, i.e.

(4.7)\begin{align} \alpha+3\beta Ohk^2&={-}\frac{3k^2(1-\beta)Oh}{1+De\alpha}- \frac{k^2E\varGamma_0}{2(1-\varGamma_0) \left(\alpha+\dfrac{k^2}{Pe}\right)}\nonumber\\ &\quad+\frac{\sigma_ik^2(1-k^2)}{2\alpha}- \frac{B_{\kappa\infty}+9B_{\mu\infty}}{2}\varGamma_0Ohk^2. \end{align}

When the Deborah number $De$ is equal to zero, (4.7) reduces to the dispersion relation determining the linear instability characteristics of surfactant-laden Newtonian viscous threads (Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018), i.e.

(4.8)\begin{equation} \alpha+3 Ohk^2={-}\frac{k^2E\varGamma_0}{2(1-\varGamma_0) \left(\alpha+\dfrac{k^2}{Pe}\right)}+\frac{\sigma_ik^2(1 -k^2)}{2\alpha}-\frac{B_{\kappa\infty}+9B_{\mu\infty}}{2}\varGamma_0Ohk^2. \end{equation}

The effect of the dimensionless numbers on the linear instability of the surfactant-covered viscoelastic thread is shown in figure 2. It can be seen that the temporal growth rate $\alpha _r$ decreases with the increase in the initial surfactant concentration $\varGamma _0$, the surface Péclet number $Pe$, the Marangoni number $E$, the dilatational Boussinesq number $B_{\kappa \infty }$ or the shear Boussinesq number $B_{\mu \infty }$ (see figures 2a–2f). In general, surfactants suppress the linear instability of viscoelastic threads, as they affect Newtonian viscous threads (Hansen et al. Reference Hansen, Peters and Meijer1999; Timmermans & Lister Reference Timmermans and Lister2002; Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018). The mechanism can be interpreted from several aspects: the presence of the surfactant lowers the surface tension and thereby the capillary force that destabilizes the viscoelastic thread; the Marangoni stress tends to eliminate the spatial non-uniformity of the surface tension and thereby suppresses the growth of the disturbance; the surface shear and dilatational viscosities and the resulting viscous stresses have a stabilizing effect on the viscoelastic thread.

Figure 2.

The temporal growth rate $\alpha _{r}$ vs the axial wavenumber $k$. The effect of (a–b) the initial surfactant concentration $\varGamma _0$, (c) the surface Péclet number $Pe$, (d) the Marangoni number $E$, (e) the dilatational Boussinesq number $B_{\kappa \infty }$, ( f) the shear Boussinesq number $B_{\mu \infty }$ and (g) the finite extensibility parameter $L$. In (a), $Pe=10^5$; in (b), $Pe=1$. (a–f) The Oldroyd-B model ($L\rightarrow \infty$) and (g) the FENE-P model. The arrows denote the direction of a parameter increasing.

Comparing the large $Pe$ case in figure 2(a) with the moderate $Pe$ case in figure 2(b), the temporal growth rate $\alpha _r$ is apparently more sensitive to the presence of the surfactant at larger values of $Pe$. Recall that the surface Péclet number $Pe$ represents the relative importance of surface convection and surface diffusion in surfactant transport. At large values of $Pe$, surface convection is dominant, resulting in high non-uniformity of surfactant and thereby large Marangoni stress that suppresses the instability of the thread. This trend can also be seen in figure 2(c). The cutoff wavenumber $k_c$, beyond which the thread is linearly stable ($\alpha _r\leq 0$), remains at unity no more in the presence of the surfactant. Analytically the cutoff wavenumber $k_c$ can be obtained from the dispersion relation (4.4). Setting the complex frequency $\alpha$ in (4.4) to be zero, we have

(4.9)\begin{equation} k_c=\sqrt{1-\frac{E\varGamma_0}{\sigma_i(1-\varGamma_0)}}. \end{equation}

Clearly, $k_c$ is strictly equal to unity only when $\varGamma _0$ is zero. As shown in figure 2, $k_c$ decreases (that is, the range of unstable axial wavenumbers is narrowed down) as $\varGamma _0, Pe$ or $E$ increases. However, this statement may not be true since the high-order terms are neglected in the 1-D model, as pointed out by Hansen et al. (Reference Hansen, Peters and Meijer1999), Timmermans & Lister (Reference Timmermans and Lister2002) and Martínez-Calvo & Sevilla (Reference Martínez-Calvo and Sevilla2018) in the study of Newtonian threads. The 2-D axisymmetric linear instability analysis has showed that the cutoff wavenumber remains at unity. The 1-D model may overestimate the suppressing effect of the surfactant on the thread instability at large wavenumbers near the cutoff. To eliminate the discrepancy, necessary high-order terms should be retained in the long wave approximation (Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018), which is of interest but beyond the scope of the present research. Anyway, the 1-D model does not need to give results that agree with the 2-D instability analysis over the entire unstable wavenumber range. We can still get some insights into how surfactants affect the linear instability and nonlinear behaviour of viscoelastic threads via the 1-D model.

A non-Newtonian viscoelastic thread possesses a higher temporal growth rate than its Newtonian viscous counterpart. That is, elasticity has a destabilizing effect on the linear instability of a thread (Brenn et al. Reference Brenn, Liu and Durst2000). It is found that finite extensibility affects viscoelastic threads in a different way. As shown in figure 2(g), the temporal growth rate $\alpha _r$ increases as the finite extensibility parameter $L$ decreases, suggesting that the extensibility of polymer chains suppresses the instability of the thread.

4.2. Nonlinear behaviour of a surfactant-covered Oldroyd-B viscoelastic thread

An Oldroyd-B viscoelastic liquid thread never breaks up because of its infinite extensibility. The 1-D models, the 2-D axisymmetric numerical simulations and the experiments have shown that a quasistatic beads-on-a-string structure is gradually formed during the nonlinear evolution of an Oldroyd-B viscoelastic thread (Chang et al. Reference Chang, Demekhin and Kalaidin1999; Ardekani et al. Reference Ardekani, Sharma and Mckinley2010; Tembely et al. Reference Tembely, Vadillo, Mackley and Soucemarianadin2012; Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018). Moreover, the filament between two adjacent primary droplets, with an almost uniform thickness, thins exponentially in time at a rate of $1/3De$, and the flow in the filament is uniaxial extensional (Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006a). The influence of other factors such as electric field, temperature, surfactant and exterior fluid phase on the dynamics of viscoelastic jets and threads has drawn attention in recent years (Eggers & Villermaux Reference Eggers and Villermaux2008; Montanero & Gañán-Calvo Reference Montanero and Gañán-Calvo2020). In the following, we examine how the presence of a surfactant monolayer affects the nonlinear behaviour of an Oldroyd-B viscoelastic thread.

The effect of the initial surfactant concentration $\varGamma _0$ on the temporal evolution of the Oldroyd-B viscoelastic thread is shown in figure 3, where the large surface Péclet number case ($Pe=10^5$) is considered and the surface dilatational and shear viscosities are turned off ($B_{\kappa \infty }=B_{\mu \infty }=0$). As shown in figure 3(a), the radius of the Oldroyd-B viscoelastic thread at the midpoint $z=\lambda /2$ (where the centre of the filament is located), denoted by $S_{\lambda /2}$, first experiences a fast decrease and then enters a slow exponential thinning, regardless of the presence of the surfactant. The surfactant has a negligible influence on the $1/3De$ exponential thinning law of the filament. A local balance between the capillary and elastic forces is presumably established in the filament as in the absence of the surfactant. The first normal stress difference $\tau _{zz}-\tau _{rr}$ increases exponentially in time at the rate $1/3De$, as shown in figure 3(b). It indicates that the surfactant plays a minimal role at the late stages of filament thinning. In contrast, the surfactant does greatly influence the evolution of the thread at the early stages. Most notably, the decrease of $S_{\lambda /2}$ is slowed down in the presence of the surfactant. As shown in figure 3(a), the filament needs increasingly more time to reach the uniaxial extension stage, as $\varGamma _0$ increases. This trend is qualitatively in accordance with the linear prediction in figure 2(a). The temporal evolution of $S_{\lambda /2}$ predicted by linear theory is included in figure 3(a) for comparison (the solid lines). The 1-D linear theory turns out to underestimate the development of the disturbance as the nonlinear effects come into play.

Figure 3.

Effect of the initial surfactant concentration $\varGamma _0$ on the thinning of the Oldroyd-B viscoelastic thread. (a) The thread radius $S$ and (b) the first normal stress difference $\tau _{zz}-\tau _{rr}$ at the midpoint $z=\lambda /2$. Here, $k=0.8$. The solid lines in (a) are the predictions of linear theory.

To understand the negligible effect of the surfactant on the thinning of the viscoelastic thread at the late stages, the variation of the thread radius $S$, the surfactant concentration $\varGamma$, the surface tension $\sigma$, the axial velocity $u$ and the first normal stress difference $\tau _{zz}-\tau _{rr}$ along the thread are examined in figure 4(a) for the large surface Péclet number case ($Pe=10^5$) and the initial surfactant concentration $\varGamma _0=0.3$ (the reference state). Obviously, at the observed time $t=212$, the beads-on-a-string structure has been formed and an elasto-capillary balance has been established inside the filament. Note that the surfactant concentration $\varGamma$ varies along the thread in a way quite similar to the thread radius $S$. As a matter of fact, at a value of the surface Péclet number as large as $10^5$, the diffusion term on the right-hand side of the surfactant advection–diffusion equation (3.7) or (3.11) is very small. Neglecting the diffusion term, (3.8) approximates to

(4.10)\begin{equation} \frac{\partial(\varGamma^2)}{\partial t}+\frac{\partial(\varGamma^2u)}{\partial z}=0.\end{equation}

The above equation is the same in form with the volume conservation equation (3.1), implying that the spatio-temporal evolution of the surfactant concentration $\varGamma$ and the thread radius $S$ follows the same discipline. Hence, we can write

(4.11)\begin{equation} \varGamma(z, t) = aS(z, t), \end{equation}

where $a$ is a coefficient to be determined. Apparently, $a=\varGamma _0$. The surface tension $\sigma$ in figure 4(a) is approximated as

(4.12)\begin{equation} \sigma= 1+E \ln(1-\varGamma_0S). \end{equation}

Most of surfactant molecules initially residing at the filament surface have been convected to the primary droplet before the thread enters the final elasto-capillary stage, due to the strong convection. The filament is almost free of surfactant in the elasto-capillary regime, with a nearly zero $\varGamma$ and a nearly unity $\sigma$. Due to the evacuation of the surfactant, the $1/3De$ exponential thinning law is maintained, the axial velocity $u$ goes linearly with $z$ and the first normal stress difference $\tau _{zz}-\tau _{rr}$ is spatially uniform in the filament, as illustrated in figure 4(a).

Figure 4.

Variation of the thread radius $S$, the surfactant concentration $\varGamma$, the surface tension $\sigma$, the axial velocity $u$ and the first normal stress difference $\tau _{zz}-\tau _{rr}$ along the thread. In panel (a), $Pe=10^5$; in panel(b), $Pe=1$. Here, $k=0.8$.

For a viscoelastic thread with no surfactant, the shape of primary droplets is almost spherical under the action of surface tension (Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006a). If the viscosity is small, primary droplets undergo shape oscillations until energy is dissipated (Ardekani et al. Reference Ardekani, Sharma and Mckinley2010; Li, Yin & Yin Reference Li, Yin and Yin2017). For the moderate viscosity considered in this problem, no oscillation takes place. On the other hand, the presence of the surfactant is found to alter the beads-on-a-string morphology to some extent. In figures 5(a) and 5(b), the thread profiles are plotted for different values of the initial surfactant concentration $\varGamma _0$ and a fixed value of the axial wavenumber $k=0.5$. It is shown that the surfactant leads to the stretch of the primary droplets in the axial direction. The primary droplets are spherical no more. They look like prolate spheroids in shape. Moreover, they deform more and more as the initial surfactant concentration $\varGamma _0$ increases. The deformation may be caused by the Marangoni stress that tends to drive liquid from the small surface tension region (the primary droplets) to the large surface tension region (the filament).

Figure 5.

(a) Thread profiles for different values of the initial surfactant concentration $\varGamma _0$, where $k=0.5$.(b) Zoomed-in secondary droplets located at the midpoint $z=\lambda /2$, where $t=100.85$ ($\varGamma _0=0$), $100.55$ ($\varGamma _0=0.1$), $105.1$ ($\varGamma _0=0.2$) and $113.9$ ($\varGamma _0=0.3$). (c) Dependence of the secondary droplet volume $V_s$ on $\varGamma _0$ and $k$.

There is a secondary droplet between two adjacent primary ones, located at the midpoint $z=\lambda /2$. Increasing the surfactant concentration results in the decrease in the secondary droplet size, as illustrated in figure 5(b). The volume of the secondary droplet, denoted by $V_s$, is shown in figure 5(c). Clearly, $V_s$ decreases as the initial surfactant concentration $\varGamma _0$ increases. When $\varGamma _0$ is large enough, the secondary droplet can be completely removed. For instance, for the axial wavenumber $k=0.7$, the secondary droplet disappears at a value of $\varGamma _0$ as large as $0.3$. The calculation result shows that the secondary droplet begins to form at the late stages of filament thinning. The mechanism, responsible for the decrease in secondary droplet size, might lie in that the tiny amount of surfactant at the filament surface slightly weakens the local capillary force and thereby impedes the formation of secondary droplet. In figure 8 below, it is shown that the Marangoni stress also plays a role in reducing the secondary droplet size. In applications, surfactants may help eliminate unwelcome secondary droplets arising in viscoelastic threads, as long as surface convection is strong enough. It is also shown in figure 5(c) that the volume of the secondary droplet increases as the axial wavenumber $k$ decreases. The most unstable wavenumber $k_{max}$ corresponding to the maximum growth rate is most probably dominant in nature. According to the 1-D linear result in figure 2(a), the most unstable wavenumber is generally small. For $\varGamma _0=0, k_{max}=0.375$; as $\varGamma _0$ increases, $k_{max}$ decreases. On the other hand, the tendency in figure 5(c) suggests that the nonlinear evolution of the viscoelastic thread exposed to natural disturbances yields secondary droplets of large volume, at least larger than those for $k=0.4$.

The moderate surface Péclet number case is examined in figure 6. At $Pe=1$, surface diffusion is comparable to surface convection and the motion of the surfactant is governed by both factors. The temporal evolution of the thread radius at the midpoint, $S_{\lambda /2}$, is illustrated in figure 6(a) for different values of the initial surfactant concentration $\varGamma _0$ and a fixed axial wavenumber $k=0.9$. For a small value of $\varGamma _0$ equal to 0.1, $S_{\lambda /2}$ decreases all the way and the thread evolves into a beads-on-a-string structure with no secondary droplet. At the large times the thinning of the filament is little influenced by the surfactant and basically follows the $1/3De$ exponential law in time. The scenario is changed with $\varGamma _0$ increasing, see the lines for $\varGamma _0$=0.3 and 0.5 in figure 6(a). At some instant, $S_{\lambda /2}$ begins to go upward, indicating that a secondary droplet is formed at the midpoint. Moreover, the size of the secondary droplet increases as $\varGamma _0$ increases. The predictions of linear theory are shown in figure 6(a) (the solid lines) for comparison.

Figure 6.

(a) Effect of the initial surfactant concentration $\varGamma _0$ on the temporal evolution of the thread radius at the midpoint, $S_{\lambda /2}$, where $Pe=1$ and the solid lines are the predictions of linear theory. Spatial distribution of (b,c) the surfactant concentration $\varGamma$ and (d,e) the surface tension $\sigma$ along the thread as time increases, where $\varGamma _0=0.3$. In (b) and (d), $Pe=10^5$; in (c) and (e), $Pe=1$. Here, $k=0.9$.

The spatio-temporal evolution of the surfactant concentration $\varGamma$ and the surface tension $\sigma$ is illustrated in figures 6(b)–6(e) for the large $Pe$ case and the moderate $Pe$ case, respectively, where the initial surfactant concentration $\varGamma _0$ is fixed to 0.3. Compared with the large $Pe$ case, in the moderate $Pe$ case the surfactant layer remains almost equally dense over the entire thread due to the surface diffusion mechanism. Recall that the surface area of the thread tends to decrease all the way to minimize the surface energy. Reasonably, the surfactant concentration $\varGamma$ increases monotonically with time from its initial value $\varGamma _0$, and the surface tension $\sigma$ decreases monotonically. Due to its abundance and non-uniform distribution at the filament surface, the surfactant affects significantly the thinning characteristics of the filament. An elasto-capillary balance fails to be established within the filament. Under the action of the capillary force and possibly the Marangoni stress as well, a small quantity of liquid is driven towards the midpoint $z=\lambda /2$ and forms a secondary droplet there.

For the case of moderate surface Péclet number $Pe=1$, the thread profiles are plotted in figure 7 for different values of the initial surfactant concentration $\varGamma _0$ and the axial wavenumber $k=0.5$. Most interestingly, the secondary droplet at the midpoint gains weight as $\varGamma _0$ increases. This trend is opposite to the large $Pe$ case shown in figure 5. In addition, high-order secondary droplets appear in the filament when $\varGamma _0$ is sufficiently large, see the profile for $\varGamma _0=0.5$.

Figure 7.

The case of moderate surface Péclet number $Pe=1$. (a) Thread profiles for different values of the initial surfactant concentration $\varGamma _0$ and (b) zoomed-in secondary droplets at the midpoint $z=\lambda /2$. Here, $k=0.5$ and $t=100.85$ ($\varGamma _0=0$), $95.75$ ($\varGamma _0=0.1$), $99.4$ ($\varGamma _0=0.3$) and $122.3$ ($\varGamma _0=0.5$).

The effect of the Marangoni stress is explored in figure 8 by comparing the case in which the Marangoni stress is taken into account with the case in which it is turned off. Figure 8(a) illustrates the temporal evolution of the thread radius at the midpoint, $S_{\lambda /2}$, for the axial wavenumbers $k=0.8$ and $0.9$. For such large axial wavenumbers, no secondary droplet is formed in the beads-on-a-string structure. It can be seen that the Marangoni stress slows down the growth of the disturbance and delays the arrival of the elasto-capillary stage to a great extent. That is, the Marangoni stress acts as a resistance to the instability of the viscoelastic thread. On the other hand, the Marangoni stress has no significant effect on the exponential thinning law of the filament in the elasto-capillary regime (recall that this is the case of large surface Péclet number). The effect of the Marangoni stress on the thread profile is examined in figures 8(b) and 8(c). For small axial wavenumbers such as $k=0.5$ and $0.6$, a secondary droplet unavoidably appears at the midpoint of the filament. The Marangoni stress decreases the size of the secondary droplet to a certain extent. In addition, the Marangoni stress stretches the primary droplets in the axial direction and makes them more prolate, as illustrated previously in figure 5.

Figure 8.

Effect of the Marangoni stress. (a) Temporal evolution of the thread radius at the midpoint, $S_{\lambda /2}$, for $k=0.8$ and $0.9$; (b) and (c) thread profiles for $k=0.5$ and $0.6$, respectively. Solid lines: the Marangoni stress is taken into account; dotted lines: the Marangoni stress is turned off. In (b), $t=113.9$ (Marangoni) and $92$ (no Marangoni); in (c), $t=126.85$ (Marangoni) and $104.6$ (no Marangoni).

The primary droplets attain an almost static shape in the elasto-capillary regime. It is of interest to seek the static solution to the shape of primary droplets. For a surfactant-free viscoelastic thread, the static solution is a sphere (Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006a), as illustrated in figure 9(a). When the thread is covered with a surfactant monolayer, its primary droplets are spherical no more. The static solution can be readily obtained in some simplified cases. For instance, if we neglect the time-dependent term, the inertia, the viscoelastic stress and the surface shear and dilatational viscosities (in the primary droplet region this simplification is reasonable), the momentum equation (3.2) yields

(4.13)\begin{equation} \frac{\partial (\sigma \mathcal{C})}{\partial z}-\frac{2}{S}\frac{\partial \sigma}{\partial z}=0. \end{equation}

In the large surface Péclet number case, the surface tension $\sigma$ can be approximated by (4.12). The substitution of (4.12) into (4.13) yields

(4.14)\begin{equation} \frac{\partial (\sigma \mathcal{C})} {\partial z}+\frac{2E\varGamma_0}{S(1-\varGamma_0S)}\frac{\partial S}{\partial z}=0. \end{equation}

Note that both $\sigma$ and $\mathcal {C}$ are functions of $S(z)$ only. So (4.14) is essentially an ordinary differential equation for $S$

(4.15)\begin{equation} \frac{\partial}{\partial z} [ \sigma \mathcal{C} +4E\varGamma_0\arctan(2\varGamma_0S-1)]=0. \end{equation}

An analytical solution to (4.14) or (4.15) may exist or not. We turn to numerical methods such as finite difference.

Figure 9.

Shape of the primary droplet located at the origin $z=0$. The cases of (a) a surfactant-free viscoelastic thread, (b) a surfactant-covered viscoelastic thread with the Marangoni stress turned off and (c) a surfactant-covered viscoelastic thread with the Marangoni stress taken into account. The dotted lines denote an isolated spherical droplet of the same volume. The solid lines denote the primary droplets simulated by the 1-D model. The dashed lines denote the static solutions. Here, $k=0.8$.

A static solution is plotted in figure 9(c). It is shown that the static solution agrees well with the 1-D numerical simulation. As outlined previously, the primary droplet exhibits a prolate spheroid profile under the action of the Marangoni stress. In contrast, when the Marangoni stress is turned off, the primary droplet becomes oblate, as illustrated in figure 9(b). That is, the primary droplet is stretched in the radial direction, due to the spatial non-uniformity of the surface tension.

The effect of the Marangoni number $E$ on the nonlinear evolution of the Oldroyd-B viscoelastic thread is shown in figure 10. As shown in figure 10(a), the Marangoni number $E$ prolongs the early stages of thread thinning. As the Marangoni number increases, the local capillary force decreases and the Marangoni stress increases. Both factors result in the decrease in the temporal growth rate of the initial disturbance, as predicted by linear theory. At the final stage, an elasto-capillary balance is established in the filament and the $1/3De$ exponential law is little influenced by the Marangoni number (recall that this is the large surface Péclet number case). On the other hand, the Marangoni number influences the beads-on-a-string morphology. Particularly, it decreases the size of the secondary droplet or even leads to the disappearance of the secondary droplet, as illustrated in figure 10(b). The capillary force at the filament surface decreases as the Marangoni number increases. This tendency might be responsible for the weight loss of the secondary droplet. The 1-D linear result in figure 2(d) shows that the most unstable wavenumber $k_{max}$ decreases from 0.375 as the Marangoni number $E$ increases from zero. Figure 10(b) shows that the smaller the axial wavenumber, the larger the secondary droplet. So secondary droplets formed under natural excitations are predicted to possess volumes larger than those for $k=0.4$ in figure 10(b).

Figure 10.

Effect of the Marangoni number $E$ on the nonlinear behaviour of the Oldroyd-B viscoelastic thread. (a) Temporal evolution of the thread radius at the midpoint, $S_{\lambda /2}$, where $k=0.8$ and the solid lines are the predictions of linear theory. (b) Dependence of the secondary droplet volume $V_s$ on $E$ and $k$. Note that for $k=0.7$ and $E=0.4$ or $0.6$ the secondary droplet actually does not exist.

The measurement of surface shear and dilatational viscosities is technically challenging due to the difficulty of separating the surface rheological effects from the Marangoni effects (Kamat et al. Reference Kamat, Wagoner, Thete and Basaran2018; Wee et al. Reference Wee, Wagoner, Kamat and Basaran2020). Several recent analytical and numerical studies on Newtonian threads have accounted for the surface viscosities and proposed possible routes to measuring the surface rheological properties (Wee et al. Reference Wee, Wagoner, Kamat and Basaran2020, Reference Wee, Wagoner, Garg, Kamat and Basaran2021). The effect of the surface viscosities on Newtonian threads are summarized in table 3. In the following, we examine their effect on the nonlinear dynamics of Oldroyd-B viscoelastic threads.

Table 3.

The effects of surface viscosities on inviscid and viscous threads.

Figure 11 illustrates the temporal evolution of the thread radius $S$, the first normal stress difference $\tau _{zz}-\tau _{rr}$, the extension rate $\dot {\varepsilon }$($=\partial u/\partial z$) and the Trouton ratio $\eta _E/\eta _0$($=3\beta +(\tau _{zz}-\tau _{rr})/Oh\dot {\varepsilon }, \eta _E$: the elongational viscosity) at the midpoint $z=\lambda /2$ for the axial wavenumber $k=0.8$, where a wide range of the Boussinesq numbers is considered. It is shown that the thinning of the viscoelastic thread can be delayed to an extraordinary extent at large values of the Boussinesq numbers, in agreement with the tendency predicted by linear theory in figures 2(e) and 2( f).

Figure 11.

Effect of the Boussinesq numbers $B_{\kappa \infty }$ and $B_{\mu \infty }$ on the thinning behaviour of the Oldroyd-B viscoelastic thread, (a) the thread radius $S$, (b) the first normal stress difference $\tau _{zz}-\tau _{rr}$, (c) the extension rate $\dot {\varepsilon }$ and (d) the Trouton ratio $\eta _E/\eta _0$ at the midpoint $z=\lambda /2$. The ordinate axis on the right is only for $B_{\kappa \infty }=B_{\mu \infty }=100$ at extremely large times. Here, $k=0.8$.

At large values of the Boussinesq numbers such as 10 and 100, the final stage of filament thinning is characterized by a uniaxial extensional flow no more. If the flow were uniaxial extensional, the extension rate $\dot {\varepsilon }$ would approach the constant $2/3De$. However, as shown in figure 11(c), for $B_{\kappa \infty }=B_{\mu \infty }=10$ or $100$, the extension rate oscillates violently with decreasing amplitude, possibly attributed to the interplay between the surface viscous stresses, the capillary force and the viscoelastic stresses as the thread enters the final stage. At the final stage, the filament thinning deviates from the $1/3De$ exponential law. Meanwhile, the temporal evolution of the first normal stress difference $\tau _{zz}-\tau _{rr}$ does not follow exactly the $1/3De$ exponential law, as shown in figure 11(b). The Trouton ratio $\eta _E/\eta _0$ behaves weirdly as well, see figure 11(d). In general, at large Boussinesq numbers, the surface viscous stresses enter the dominant force balance and affect the thinning of the filament significantly.

An interesting phenomenon is observed at moderate Boussinesq numbers. As shown in figure 11(a), for $B_{\kappa \infty }=B_{\mu \infty }=1$, the thread radius at the midpoint, $S_{\lambda /2}$, begins to increase at some instant, indicating that a secondary droplet is formed. This is a delicate situation, for the secondary droplet disappears when the Boussinesq numbers increase a little. There seems to be a narrow window of $B_{\kappa \infty }$ and $B_{\mu \infty }$ near $O(1)$ for the appearance of secondary droplet. The surface viscosities prompt the formation of secondary droplet when they are small but suppress it when they are large.

The role of the surface viscous stresses at the final stage of thread thinning can be better understood by defining the filament Ohnesorge number $Oh^S_{1f}=\mu _s(\rho \sigma R^3_f)^{-1/2}$, where the surface shear viscosity $\mu _s$, the surface tension $\sigma$ and the filament radius $R_f$ are all evaluated at the midpoint $z=\lambda /2$. The filament Ohnesorge number $Oh^S_{1f}$ represents the ratio of the surface viscous stresses to the capillary force in the filament. The evolution of $Oh^S_{1f}$ in time is shown in figures 12(a) and 12(b) for $B_{\kappa \infty }=B_{\mu \infty }=0.1$ and $B_{\kappa \infty }=B_{\mu \infty }=10$, respectively, where the thread radius $S_{\lambda /2}$ and the surfactant concentration $\varGamma _{\lambda /2}$ are plotted as well. As illustrated in the figure, before the final stage, $Oh^S_{1f}$ increases slowly and $S_{\lambda /2}$ and $\varGamma _{\lambda /2}$ decrease rapidly with time increasing; in contrast, at the final stage, $Oh^S_{1f}$ increases rapidly and $S_{\lambda /2}$ and $\varGamma _{\lambda /2}$ decrease slowly with time increasing. For the large surface Péclet number case ($Pe=10^5$) considered in figure 12, the increase in $Oh^S_{1f}$ at the final stage, similar to the decrease in $S_{\lambda /2}$ or $\varGamma _{\lambda /2}$, seems to follow an almost exponential law in time, especially for small Boussinesq numbers. For $B_{\kappa \infty }=B_{\mu \infty }=0.1, Oh^S_{1f}$ can reach $O$(1). For $B_{\kappa \infty }=B_{\mu \infty }=10, Oh^S_{1f}$, which is already larger than unity at the initial time, can be $O$($10^2$) at the final time. The surface viscous stresses can still be large in the filament and compete with the capillary force and the viscoelastic stresses, even though surfactant molecules are almost swept away from the filament for large surface Péclet numbers.

Figure 12.

Temporal evolution of the filament Ohnesorge number $Oh^S_{1f}$, the thread radius $S$ and the surfactant concentration $\varGamma$ at the midpoint $\lambda /2$. In (a), $B_{\kappa \infty }=B_{\mu \infty }=0.1$; in (b), $B_{\kappa \infty }=B_{\mu \infty }=10$. Here, $k=0.8$.

The effect of the surface viscosities on the morphology of the Oldroyd-B viscoelastic thread is further examined in figure 13, where the thread profiles are plotted for different values of the axial wavenumber $k$. The two distinct roles of the surface viscosities can be more clearly seen in the figure. For instance, for $k=0.7$, there exists a secondary droplet at the midpoint when the surface viscosities are absent ($B_{\kappa \infty }=B_{\mu \infty }=0$).With the Boussinesq numbers $B_{\kappa \infty }$ and $B_{\mu \infty }$ increasing from zero, the secondary droplet is continuously shrunk. When the surface viscosities are sufficiently large, the secondary droplet completely disappears, see the thread profile for $B_{\kappa \infty }=B_{\mu \infty }=100$, reflecting the suppression effect of the surface viscosities on secondary droplets. The same phenomenon occurs at smaller axial wavenumbers, see the thread profiles for $k=0.3$ to $0.6$. This trend is easily understood, considering that large surface viscosities create large surface viscous stresses that tend to hinder the formation of secondary droplets. For larger axial wavenumbers $k=0.8$ and $0.9$, the thread evolves into a beads-on-a-string structure with no secondary droplet, in the absence of surface viscosities. However, at $B_{\kappa \infty }=B_{\mu \infty }=1$, the surface viscosities cause the appearance of a secondary droplet at the midpoint. Then, the surface viscosities try to shrink the secondary droplet and eventually lead to its disappearance, see the thread profiles for $B_{\kappa \infty }=B_{\mu \infty }=100$. The formation of secondary droplet at large axial wavenumbers and moderate surface viscosities seems odd. It possibly results from the interplay between the inertial, capillary and surface viscous forces.

Figure 13.

Thread profiles for different values of the Boussinesq numbers $B_{\kappa \infty }$ and $B_{\mu \infty }$ and different values of the axial wavenumber $k$. The insets illustrate the zoomed-in secondary droplet located at the midpoint ${z=\lambda /2}$.

It is also found that surface viscosities can deform primary droplets. As shown in figure 13, the primary droplets are continuously stretched in the axial direction as the shear and dilatational Boussinesq numbers increase. At a value of $B_{\kappa \infty }$ and $B_{\mu \infty }$ as large as $100$, the shape of the primary droplets becomes singular. The corner regions connecting the filament to the primary droplets get stretched. The shape resembles a cone in the 3-D space. From a physical point of view, the capillary force fails to overcome the surface viscous stresses and cannot help the primary droplets maintain their spheroidal shape at high surface viscosities.

In the 1-D momentum equation (3.2), the terms numbered ${\bigcirc{\kern-6pt 2}}$, ${\bigcirc{\kern-6pt 4}}$ and ${\bigcirc{\kern-6pt 5}}$ represent the surface viscous stresses. If the shear and dilatational viscosities are not equal, terms ${\bigcirc{\kern-6pt 2}}$ and ${\bigcirc{\kern-6pt 4}}$ are non-zero. These two terms make the simulation more difficult. To simplify the calculation, one strategy is to apply the slender body approximation to the full curvature in ${\bigcirc{\kern-6pt 2}}$ and ${\bigcirc{\kern-6pt 4}}$, i.e.

(4.16)\begin{equation} \mathcal{C}\simeq \frac{1}{S}. \end{equation}

The accuracy of this simplification is checked in figure 14 by comparing the slender body approximation with the full-curvature case. It is shown in figure 14(a) that the simplification does not lead to any significant difference in the temporal evolution of the Oldroyd-B viscoelastic thread. In figures 14(b) and 14(c), the slender body approximation seems to predict smaller secondary droplets and less spherical primary droplets. Note that $B_{\kappa \infty }$ is taken to be smaller than $B_{\mu \infty }$. If $B_{\kappa \infty }$ is larger than $B_{\mu \infty }$, the trends in figures 14(b) and 14(c) are presumably opposite, since the signs of $B_{\kappa \infty }$ and $B_{\mu \infty }$ in terms ${\bigcirc{\kern-6pt 2}}$ and ${\bigcirc{\kern-6pt 4}}$ of (3.2) are opposite.

Figure 14.

Comparison of the slender body approximation (solid lines) and the full-curvature (dotted lines) case. (a) Temporal evolution of the thread radius at the midpoint, $S_{\lambda /2}$, for $k=0.8$ and $0.9$, where $B_{\kappa \infty }=1$ and $B_{\mu \infty }=10$; (b) and (c) thread profiles for $k=0.7$ and $k=0.8$, respectively, where $B_{\kappa \infty }=0.1$ and $B_{\mu \infty }=1$.

4.3. Nonlinear behaviour of a surfactant-covered FENE-P viscoelastic thread

Different from Oldroyd-B viscoelastic threads, viscoelastic threads of FENE-P type break up in finite time, owing to finite extensibility of polymer chains. The nonlinear behaviour of a FENE-P viscoelastic thread free of surfactant is represented in figure 15. As shown in figure 15(a), the $1/3De$ exponential thinning law is valid no more for the FENE-P model. At the nonlinear stages, the elastic force in the filament fails to balance the capillary force, and the dynamics of the thread is dominated by inertia, capillarity and viscoelasticity. In such a case, the minimum radius of the thread decreases linearly in time until pinch-off. For smaller values of the finite extensibility $L$, the minimum thread radius $S_{min}$ decreases faster and pinch-off occurs earlier. Meanwhile, the location of $S_{min}$, denoted by $z_{min}$, moves from the midpoint $z_{\lambda /2}$ to the corner region connecting the filament to the primary droplet, although the thickness of the filament appears to be uniform (Fontelos & Li Reference Fontelos and Li2004). The pinch-off time, denoted by $t_p$, is diagrammed in figure 15(b). It can be seen that $t_p$ decreases as $L$ decreases. On the other hand, $t_p$ depends greatly on the axial wavenumber, or rather, on the growth rate predicted by linear theory. The thread profiles are plotted in figures 15(c) and 15(d) for a relatively large axial wavenumber $k=0.8$ and a relatively small axial wavenumber $k=0.3$, respectively. It is found that finite extensibility may alter the beads-on-a-string morphology to a certain extent. In the case of $k=0.8$, a tiny secondary droplet arises at the midpoint when $L$ is sufficiently small, say $L=2$. In the case of $k=0.3$, small extensibility may lead to large deformation of secondary droplet, see the thread profile for $L=5$. In the following, we examine briefly the effect of surfactants on the nonlinear behaviour of FENE-P viscoelastic liquid threads.

Figure 15.

Effect of the finite extensibility parameter $L$ on the nonlinear behaviour of a FENE-P viscoelastic thread free of surfactant. (a) The minimum thread radius $S_{min}$ (the lower lines) and its location $z_{min}$ (the upper lines) for different values of $L$, where $k=0.8$. (b) The pinch-off time $t_p$ vs $L$ for different values of $k$. Thread profiles for (c) $k=0.8$ and (d) $k=0.3$ for different values of $L$.

Figure 16 illustrates the effect of the initial surfactant concentration $\varGamma _0$ on the thinning of the FENE-P viscoelastic thread. It can be seen from figure 16(a) that the surfactant retards the pinch-off of the thread. This is mainly because the decrease in the minimum thread radius $S_{min}$ at the early stages is decelerated with increasing $\varGamma _0$, which is basically a linear effect. At the final stage, the influence of the surfactant is limited, for most of the surfactant molecules are swept away from the filament for the large surface Péclet number case considered, as discussed in the Oldroyd-B model. In figure 16(b), close to the pinch-off time $t_p, S_{min}$ goes to zero almost in the same way as $\varGamma _0$ varies, only that $S_{min}$ decreases a little faster at larger values of $\varGamma _0$.

Figure 16.

Effect of the initial surfactant concentration $\varGamma _0$ on the thinning of the FENE-P viscoelastic thread. (a) Temporal evolution of the minimum thread radius $S_{min}$ and (b) $S_{min}$ vs the time to the pinching $t_p-t$ for different values of $\varGamma _0$. Here, $k=0.8$.

The effect of surfactant surface diffusion on the nonlinear behaviour of the FENE-P viscoelastic thread is explored by considering the moderate surface Péclet number case $Pe=1$. The result is illustrated in figure 17. As shown in figure 17(a), similar to the large $Pe$ case ($Pe=10^5$), in the moderate $Pe$ case the pinch-off of the FENE-P viscoelastic thread is retarded by increasing the initial surfactant concentration $\varGamma _0$, owing to the suppression effect of the surfactant on the temporal growth rate of the initial disturbance. At the nonlinear stages, the minimum thread radius $S_{min}$ approaches zero at different speeds for different values of $\varGamma _0$, as illustrated in figure 17(b). Unfortunately, no common trend can be found, suggesting that the surfactant affects the pinch-off of the viscoelastic thread in a complicated way.

Figure 17.

The case of moderate surface Péclet number $Pe=1$. (a) Temporal evolution of the minimum thread radius $S_{min}$ and (b) $S_{min}$ vs the time to the pinching $t_p-t$ for distinct values of the initial surfactant concentration $\varGamma _0$. (c) Spatial variation of the thread radius $S$, the surfactant concentration $\varGamma$, the surface tension $\sigma$, the axial velocity $u$ and the first normal stress difference $\tau _{zz}-\tau _{rr}$ along the thread. In (a–c), $k=0.8$. In (c), the solid lines are for $\varGamma _0=0.3$ (the time instant $t=110.795$) and the dotted lines are for $\varGamma _0=0$ ($t=94.48$).(d) Thread profiles for $k=0.8$ and distinct values of $\varGamma _0$, where the inset is the zoomed-in filament and secondary droplet regions. (e) Thread profiles for $k=0.4$ and distinct values of $\varGamma _0$.

The variation of the thread radius $S$, the surfactant concentration $\varGamma$, the surface tension $\sigma$, the axial velocity $u$ and the first normal stress difference $\tau _{zz}-\tau _{rr}$ along the thread is illustrated in figure 17(c) for the axial wavenumber $k=0.8$ and the initial surfactant concentration $\varGamma _0=0$ and $0.3$. The spatial variation of $u$ and $\tau _{zz}-\tau _{rr}$ suggests that a tiny secondary droplet arises at the midpoint $z=\lambda /2$ in the absence of the surfactant ($\varGamma _0=0$), although it is almost invisible from the thread profile. In the presence of the surfactant, i.e. $\varGamma _0=0.3$, the volume of the secondary droplet is visibly increased. Due to surface diffusion, surfactant molecules cover the entire thread densely. Particularly in the neighbourhood of the secondary droplet, the surfactant concentration $\varGamma$ is close to unity and the surface tension $\sigma$ is very small. Correspondingly, the viscoelastic force in the neighbourhood of the secondary droplet is small and the secondary droplet deforms greatly from spherical.

The effect of the initial surfactant concentration $\varGamma _0$ on the beads-on-a-string morphology in the moderate $Pe$ case can be seen from figures 17(d) and 17(e) for the relatively large axial wavenumber $k=0.8$ and the relatively small axial wavenumber $k=0.4$, respectively. For both wavenumbers, the secondary droplet at the midpoint gains weight and exhibits a more singular shape as $\varGamma _0$ increases. Meanwhile, the thickness of the filament becomes more non-uniform. Pinch-off always occurs at the corner region where the filament joins the primary droplet. For the relatively small axial wavenumber $k=0.4$, higher-order generations of secondary droplets can be formed at sufficiently large values of $\varGamma _0$ (see the thread profile for $\varGamma _0=0.5$).

The effect of the Marangoni stress on the pinch-off of the FENE-P viscoelastic thread is examined in figure 18, by comparing the case in which the Marangoni stress is taken into account with the case in which the Marangoni stress is turned off. As seen in figure 18(a), for the relatively large axial wavenumber $k=0.8$, the pinch-off of the viscoelastic thread occurs much earlier when the Marangoni stress is turned off. Note that the two lines intersect with each other near $S_{min}=0.53$ and $t_p-t=23$. Before intersecting, the minimum thread radius $S_{min}$ decreases much more rapidly in the absence of the Marangoni stress. In contrast, from the intersection point, it takes $S_{min}$ the same time in two cases to decrease to zero. In general, the Marangoni stress greatly delays the pinch-off of the FENE-P viscoelastic thread. The thread profiles prior to pinch-off are plotted in figures 18(b) and 18(c) for $k=0.8$ and $0.3$, respectively. For both axial wavenumbers, the primary droplets have been axially stretched by the Marangoni stress before the pinch-off takes place. The secondary droplets remain singular in shape when the Marangoni stress is absent. It indicates that the Marangoni stress has little to do with the non-spherical shape of secondary droplets. The deformation of secondary droplets possibly results from small local surface tension. On the other hand, the Marangoni stress affects secondary droplets by reducing their size.

Figure 18.

Effect of the Marangoni stress on the nonlinear behaviour of the FENE-P viscoelastic thread.(a) The minimum thread radius $S_{min}$ vs the time to the pinching $t_p-t$ for $k=0.8$. Thread profiles for(b) $k=0.8$ and (c) $k=0.3$.

Figure 19 shows the effect of the Marangoni number $E$ on the pinch-off of the FENE-P viscoelastic thread. At the early stages, the thinning of the thread can be greatly slowed down by increasing the value of $E$, owing to the suppression effect of the Marangoni number on the linear instability of the thread. At the late stages, from approximately $S_{min}=0.45$ and $t_p-t=21$, the effect of the Marangoni number becomes quite limited. Recall that in the large surface Péclet number case of the Oldroyd-B viscoelastic thread, surfactant molecules are continuously convected to the primary droplet and not many remain at the filament surface. The situation is similar for the FENE-P viscoelastic thread.

Figure 19.

Effect of the Marangoni number $E$ on the pinch-off of the FENE-P viscoelastic thread. The minimum thread radius $S_{min}$ vs the time to the pinching $t_p-t$. Here, $k=0.8$.

In addition to the Marangoni number $E$, the Boussinesq numbers $B_{\kappa \infty }$ and $B_{\mu \infty }$ may greatly delay the pinch-off of the FENE-P viscoelastic thread, as illustrated in figure 20(a). At relatively small values of $B_{\kappa \infty }$ and $B_{\mu \infty }$ such as 0.1 and 1, the surface viscosities slow down the growth of the initial disturbance at the early stages. At relatively large values of $B_{\kappa \infty }$ and $B_{\mu \infty }$ 10 and 100, the surface viscosities tend to inhibit the pinch-off of the viscoelastic thread. The thread profiles for the axial wavenumber $k=0.3$ and several values of $B_{\kappa \infty }$ and $B_{\mu \infty }$ are plotted in figure 20(b). Like in Oldroyd-B threads, the surface viscous stresses may lead to the axial stretch of primary droplets. When $B_{\kappa \infty }$ and $B_{\mu \infty }$ are large, e.g. $B_{\kappa \infty }=B_{\mu \infty }=10$, the corner connecting the filament to the primary droplet is axially stretched and the primary droplet is greatly deformed. On the other hand, for this small axial wavenumber, the secondary droplet located at the midpoint loses weight under the action of the surface viscous stresses. At sufficiently large values of $B_{\kappa \infty }$ and $B_{\mu \infty }$, the secondary droplet disappears.

Figure 20.

Effect of the Boussinesq numbers $B_{\kappa \infty }$ and $B_{\mu \infty }$ on the pinch-off of the FENE-P viscoelastic thread. (a) The minimum thread radius $S_{min}$ vs the time to the pinching $t_p-t$ for $k=0.8$. (b) Thread profiles for different values of $B_{\kappa \infty }$ and $B_{\mu \infty }, k=0.3$.