2.1. Governing equations

First, the electrostatic effects are governed by the Maxwell equations in the absence of magnetic fields,

(2.1) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }(\varepsilon _i\boldsymbol{E}_i)=\rho _i^e,\quad \boldsymbol{\nabla }\times \boldsymbol{E}_i=0,\quad i=1,2, \end{equation}

where $\rho _i^e$ is the volume charge density within the fluid, and $\boldsymbol{E}_i$ is the electric field. We then introduce the electric potential $\phi _i$ as follows:

(2.2) \begin{equation} \boldsymbol{E}_i=-\boldsymbol{\nabla }\phi _i. \end{equation}

At the interface, electric potential must satisfy the following jump condition:

(2.3) \begin{equation} \left \| \varepsilon \boldsymbol{\nabla }\phi _i\boldsymbol{\cdot }\boldsymbol{n} \right \|=-q,\quad \left \| \boldsymbol{\nabla }\phi _i\boldsymbol{\cdot }\boldsymbol{t} \right \|=0, \end{equation}

where $q$ represents the surface charge density. The vectors $\boldsymbol{n}$ and $\boldsymbol{t}$ denote the unit normal (pointing outward) and unit tangential vectors, respectively. The notation $\left \| (\boldsymbol{\cdot }) \right \|$ represents the jump of the physical quantity $(\boldsymbol{\cdot })$ across the interface, defined as $(\boldsymbol{\cdot })_2-(\boldsymbol{\cdot })_1$ .

The electrolyte used to induce the ionic surfactant is a symmetric $z:z$ electrolyte. We use superscripts ` $+$ ’ and ` $-$ ’ to denote the positive and negative ions in the core flow, respectively. The concentration of positive and negative ions coming from the dissolution of surfactant in the liquid thread are represented by $c^{(+)}$ and $c^{(-)}$ , and the positive charge ions are adsorped into the interface and act as the surfactant, while $n^{(+)}$ and $n^{(-)}$ denote the ions of other electrolytes unaffected by the surfactant. The units of $c^{(\pm )}$ and $n^{(\pm )}$ are in molar concentration. The volume charge density in the core flow can then be expressed as

(2.4) \begin{equation} \rho _1^e={\textit{eN}}_A\big(n^{(+)}+c^{(+)}-n^{(-)}-c^{(-)}\big), \end{equation}

where $e$ is the charge on an electron, and $N_A$ is Avogadro’s number. And it is assumed that there is no charge in the air, i.e. $\rho ^e_2=0$ .

The surface charge density is given by

(2.5) \begin{equation} q={\textit{eN}}_A\big(\varGamma +n_s^{(+)}-n_s^{(-)}\big), \end{equation}

where $\varGamma$ represents the positively charged ionic surfactant concentration at the surface. And the surface species $n^{+}_s, n^{-}_s$ from other electrolytes satisfy the Nernst–Plank equation at the interface and their maximum electrolyte packing concentration is assumed to be $n^0_s$ (see (8), (9) and related discussion in Conroy et al. (Reference Conroy, Matar, Craster and Papageorgiou2011)).

The ion concentration in the core flow evolves according to the Nernst–Planck equation,

(2.6) \begin{equation} \frac {{\rm D}c^{(\pm )}}{{\rm D}t}=\boldsymbol{\nabla }\boldsymbol{\cdot }{\omega \big[ec^{(\pm )}\boldsymbol{\nabla }\phi _i+k_BT\boldsymbol{\nabla }c^{(\pm )}\big]}, \quad \frac {{\rm D}n^{(\pm )}}{{\rm D}t}=\boldsymbol{\nabla }\boldsymbol{\cdot }{\omega \big[e n^{(\pm )}\boldsymbol{\nabla }\phi _i+k_BT\boldsymbol{\nabla }n^{(\pm )}\big]}, \end{equation}

where $k_B$ is the Boltzmann constant, $T$ is the absolute temperature, $\omega$ represents the mobility of the species, ${\rm D}/{\rm D}t$ is the material derivative. Due to considering surfactants and other electrolytes as strong electrolytes, there are no dissociation and association terms for ions in (2.6). In addition, the equilibrium condition in the thread satisfies where the electric potential $\phi _1=\phi ^0$ (assuming $\phi ^0=0$ ) and the corresponding total positive and negative ion concentrations are both $c_0$ .

At the interface, the concentration $\varGamma$ of the positively charged surfactants evolves according to the following evolution equation:

(2.7) \begin{equation} \frac {\partial \varGamma }{\partial t}-\frac {\partial \boldsymbol{X}}{\partial t}\boldsymbol{\nabla\!}_s\varGamma +\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }(\varGamma \boldsymbol{u}_t)+\varGamma \zeta \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }(D_s\boldsymbol{\nabla\!}_s\varGamma +D_{s\phi }\varGamma \boldsymbol{\nabla\!}_s\phi _s), \end{equation}

where $\partial \varGamma /\partial t$ is the time derivative following the evolution of the surface point by keeping the surface coordinates fixed, $\boldsymbol{X}$ is the position vector of the surface point, $\boldsymbol{\nabla\!}_s=(\boldsymbol{I}-\boldsymbol{nn})\boldsymbol{\cdot }\boldsymbol{\nabla}$ is the surface gradient operator, $\boldsymbol{u}$ is the velocity of the surface point and $\boldsymbol{u}_t$ is the tangential component of velocity for the corresponding surface point, $\zeta =\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{n}$ is the curvature at the interface, $D_s$ is the surface diffusion coefficient, $D_{s\phi }=D_s e/k_B T$ is the diffusion coefficient for ion migration and $\phi _s$ is the electric potential at the interface. The above evolution equation is based on the book chapter (Papageorgiu Reference Papageorgiu2002), where the ion electromigration term is obtained by assuming $J_{\phi }\,D_{s\phi }\boldsymbol{\nabla\!}_s\phi _s$ . And here we make the equilibrium assumption for low concentrations and ignore surface flux to the interface from the bulk fluid (Conroy et al. Reference Conroy, Matar, Craster and Papageorgiou2011). This study considers the case where the surfactant concentration does not reach the critical micelle concentration.

The continuity equation and Navier–Stokes equations for the fluid can be written as follows:

(2.8) \begin{align}&\qquad\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_i=0,\quad i=1,2, \end{align}

(2.9) \begin{align}& \rho _i\frac {{\rm D}\boldsymbol{u}_i}{{\rm D}t}=\boldsymbol{\nabla }\boldsymbol{\cdot }\big(\boldsymbol{T}_i^f+\boldsymbol{T}_i^e\big),\quad i=1,2. \end{align}

Assuming the fluid is incompressible, the fluid stress tensor and Maxwell stress tensor appearing in (2.9) can be given as follows:

(2.10) \begin{align}&\qquad \boldsymbol{T}_i^f=-p_i\boldsymbol{I}+\boldsymbol{\tau }_i, \end{align}

(2.11) \begin{align}& \boldsymbol{T}_i^e=\varepsilon _i\left(\boldsymbol{E}_i\boldsymbol{E}_i-\frac {1}{2}|\boldsymbol{E}_i|^2\boldsymbol{I}\right)\!. \end{align}

Here, $p_1$ is pressure of the liquid thread, $p_2$ is air pressure, $\boldsymbol{\tau }_i=\mu _i[\boldsymbol{\nabla }\boldsymbol{u}_i+(\boldsymbol{u}_i)^{\textrm {{T}}}]$ and $\boldsymbol{\nabla }\boldsymbol{\cdot }T_i^e=\rho _i^e\boldsymbol{E}_i$ .

The boundary conditions can be summarized as kinematic boundary conditions and dynamic boundary conditions. The first is the no-slip boundary condition, with $\boldsymbol{u}_2(R_2,z,t)=\textbf {0}$ at the pipe wall. In addition, there is a continuity of velocity at the gas–liquid interface,

(2.12) \begin{equation} \left \| \boldsymbol{u}_i \right \|=0, \end{equation}

where $\boldsymbol{u}_i=(u_i, w_i)^T$ , $w_i$ represents the axial velocity ( $z$ direction) and $u_i$ represents the radial velocity ( $r$ direction).

The movement of the position of the gas–liquid interface conforms to the kinematic boundary conditions,

(2.13) \begin{equation} u_i=S_t+w_iS_z, \end{equation}

where the velocity term is calculated at the interface. $S_t$ denotes the partial derivative of $S$ with respect to time and $S_z$ denotes the partial derivative of $S$ with respect to $z$ .

In this model, it is assumed that surfactants induce both surface elasticity and surface viscosity at the gas–liquid interface. In other words, it not only reduces surface tension but also causes additional viscous stresses at the interface, as described by Martínez-Calvo & Sevilla (Reference Martínez-Calvo and Sevilla2018). The free-surface stress is simulated by two surface viscosities: surface shear viscosity $\mu _s(\varGamma )$ and dilatational viscosity $\kappa _s(\varGamma )$ . Following the derivation process of Martinez-Calvo & Sevilla (Reference Martínez-Calvo and Sevilla2018), the stress balance at the gas–liquid interface can be expressed as follows:

(2.14) \begin{equation} \big \| \boldsymbol{T}_i^f+\boldsymbol{T}_i^e \big \|\boldsymbol{\cdot }\boldsymbol{n}+\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }\boldsymbol{T\!}_s=0, \end{equation}

where the surface stress tensor $\boldsymbol{T\!}_s$ is defined using the well-known Boussinesq–Scriven approximation, which neglects surface density and simplifies the complex interfacial rheology. This approximation assumes that the deviatoric components of $\boldsymbol{T\!}_s$ are isotropic, linear and instantaneous in the surface strain rate tensor (Boussinesq Reference Boussinesq1913; Scriven & Sternling Reference Scriven and Sternling1960),

(2.15) \begin{equation} \boldsymbol{T\!}_s=\boldsymbol{I}_s[\sigma +(\kappa _s-\mu _s)(\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }\boldsymbol{u}_s)]+\mu _s\left[(\boldsymbol{\nabla\!}_s\boldsymbol{u}_s)\boldsymbol{\cdot }\boldsymbol{I}_s+\boldsymbol{I}_s(\boldsymbol{\nabla\!}_s\boldsymbol{u}_s)^{\textrm {T}}\right]\!. \end{equation}

Substituting the (2.15) into (2.14) yields

(2.16) \begin{equation} \begin{aligned} &\big \| \boldsymbol{T}_i^f+\boldsymbol{T}_i^e \big \|\boldsymbol{\cdot }\boldsymbol{n}+\boldsymbol{\nabla\!}_s\sigma -\boldsymbol{n}(\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }\boldsymbol{n})\sigma +\boldsymbol{\nabla\!}_s[(\kappa _s-\mu _s)(\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }\boldsymbol{u}_s)] \\ &\quad -\boldsymbol{n}(\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }\boldsymbol{n})(\kappa _s-\mu _s)(\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }\boldsymbol{u}_s)+\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }\left \{ {\mu _s\left[(\boldsymbol{\nabla\!}_s\boldsymbol{u}_s)\boldsymbol{\cdot }\boldsymbol{I}_s+\boldsymbol{I}_s\boldsymbol{\cdot }(\boldsymbol{\nabla\!}_s\boldsymbol{u}_s)^{\textrm {T}}\right]} \right \}=0, \end{aligned} \end{equation}

where $\zeta =\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }\boldsymbol{n}= {1}/{S\sqrt {1+S_z^2}}- {S_{\textit{zz}}}/{\sqrt {(1+S_z^2)^3}}$ is mean curvature of the interface.

The velocity field on the surface of the fluid, $\boldsymbol{u}_s$ , can be decomposed into two components: tangential velocity $\boldsymbol{u}_t$ and normal velocity $\boldsymbol{u}_n$ ,

(2.17) \begin{equation} \boldsymbol{u}_s=\boldsymbol{u}(\boldsymbol{X},t)=\boldsymbol{u}_t(\boldsymbol{X},t)+\boldsymbol{u}_n(\boldsymbol{X},t)=\boldsymbol{u}_t(\boldsymbol{X},t)+u_n(\boldsymbol{X},t)\boldsymbol{n}, \end{equation}

where $\boldsymbol{u}_t=\boldsymbol{I}_s\boldsymbol{\cdot }\boldsymbol{u}_s$ , $\boldsymbol{u}_n=\boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{u}_s$ and $\boldsymbol{I}_s=\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n}$ is surface projection operator. Since the interface is axisymmetric, we have $\boldsymbol{u}_t=u_t(\boldsymbol{X},t)\boldsymbol{t}$ , and then the normal and tangential velocity at the interface can be written as follows:

(2.18) \begin{equation} u_n=\frac {u-S_z w}{\sqrt {1+S_z^2}},\quad u_t=\frac {w+S_z u}{\sqrt {1+S_z^2}}, \end{equation}

where $\boldsymbol{u}_s=(u,w)^{\text{T}}$ , $u$ and $w$ are the radial and axial components of fluid velocity on the interface, respectively. In light of the aforementioned equations, $\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }\boldsymbol{u}_s=\boldsymbol{\nabla\!}_s\boldsymbol{\cdot }\boldsymbol{u}_t+\zeta u_n$ can be applied to (2.16).

Taking the inner product for both sides of (2.16) with $\boldsymbol{n}$ and $\boldsymbol{t}$ , respectively, yields the normal and tangential boundary conditions at the interface,

(2.19) \begin{equation} \begin{aligned} &-\!\left \| p_i-\frac {2\mu _i}{1\!+\!S_z^2}\left[u_r-S_z(u_z+w_r)+\!S_z^2w_z\right] \right \|\!+\!\left \| \frac {\varepsilon _i}{2}\frac {1\!-\!S_z^2}{1\!+\!S_z^2}(\phi _{\textit{ir}}^2-\phi _{\textit{iz}}^2)-2\varepsilon _i\frac {S_z^2}{1\!+\!S_z^2}\phi _{\textit{ir}}\phi _{\textit{iz}} \right \|\\&=\zeta \!\left [\! \sigma +(\kappa _s-\mu _s)\!\left ( \!\frac {({\textit{Su}}_t)_z}{S\sqrt {1\!+\!S_z^2}}+\zeta u_n \!\right )\! \right ] +\frac {2\mu _s}{1\!+\!S_z^2}\!\left [ \!\frac {S_z u_t\!+u_n}{S^2}-\frac {S_{\textit{zz}}}{1\!+\!S_z^2}\!\left ( \!\frac {\partial u_t}{\partial z}-\frac {S_{\textit{zz}}u_n}{1\!+\!S_z^2} \!\right ) \!\right ]\!, \end{aligned} \end{equation}

(2.20) \begin{align} &-\frac {1}{\sqrt {1+S_z^2}}\left \| 2\mu _i\left [ u_r S_z+\frac {1}{2}(1-S_z^2)(u_z+w_r)-S_z w_z \right ] \right \|+\frac {2\beta }{\sqrt {1+S_z^2}} q(S_z\phi _{i,r}+\phi _{i,z})|_{S} \nonumber\\& =\frac {\partial \sigma }{\partial z}+\frac {\partial }{\partial z}\!\left [\! (\kappa _s-\mu _s)\!\left ( \!\frac {({\textit{Su}}_t)_z}{S\sqrt {1\!+\!S_z^2}}+\zeta u_n \!\right )\! \right ]+\frac {\sqrt {1\!+\!S_z^2}}{S}\frac {\partial }{\partial z}\!\left [ \!\frac {2\mu _s S}{\sqrt {1\!+\!S_z^2}}\!\left ( \!\frac {\partial u_t}{\partial z}-\frac {S_{\textit{zz}}u_n}{1\!+\!S_z^2} \!\right ) \!\right ] \nonumber\\ &\quad -\frac {2\mu _s S_z}{\sqrt {1+S_z^2}}\left [ \frac {S_z u_t+u_n}{S^2}-\frac {S_{\textit{zz}}}{1+S_z^2}\left ( \frac {\partial u_t}{\partial z}-\frac {S_{\textit{zz}}u_n}{1+S_z^2} \right )\right ]\!, \\[9pt] \nonumber \end{align}

where the subscripts ‘ $r$ ’ and ‘ $z$ ’ represent the partial derivative.

It should be noted that both the surface shear viscosity $\mu _s(\varGamma )$ and the dilatational viscosity $\kappa _s(\varGamma )$ are functions of the surfactant concentration $\varGamma$ at the interface. It can be assumed that $\mu _s(\varGamma )$ and $\kappa _s(\varGamma )$ are linearly dependent on $\varGamma$ (Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018), that is

(2.21) \begin{equation} \kappa _s(\varGamma )=\kappa _{\textit{sr}}\frac {\varGamma }{\varGamma _r},\quad \mu _s(\varGamma )=\mu _{\textit{sr}}\frac {\varGamma }{\varGamma _r}, \end{equation}

where the $\varGamma _r$ represents the reference point for the surfactant concentration, $\kappa _{\textit{sr}}$ and $\mu _{\textit{sr}}$ are two constants.

Furthermore, expression of surface tension under an electrical field is elucidated, which is a function of the surfactant concentration and the interfacial electric potential. Prosser & Franses (Reference Prosser and Franses2001) has provided the formula for calculating the surface tension in the case of over equilibrium, utilizing the chemical potential and the Gibbs adsorption equation. Datwani & Stebe (Reference Datwani and Stebe1999) has also delineated the impact of an excess Maxwell stress on the surface tension generated by the electric potential on the surface. The aforementioned studies have demonstrated that the surface tension can be expressed as the sum of the contributions from the Langmuir isotherms and the electrical effect,

(2.22) \begin{equation} \sigma =\sigma _0+R_g T\varGamma _{\infty }\rm {ln}\left(1-\frac {\varGamma }{\varGamma _{\infty }}\right)+\int _0^{\phi _s}qd\phi , \end{equation}

where $\sigma _0$ denotes the surface tension of the cleaning interface, $\varGamma _{\infty }$ denotes the maximum packing surfactant concentration and $R_g$ is the gas constant.

2.2. One-dimensional model

In light of the elongated thread, the radial reference size $R_0$ is considerably smaller than the axial reference size $L$ , rendering the one-dimensional model analysis a viable approach. As the outside of the thread is air, whose viscosity is much smaller than the liquid, and air pressure can be assumed to be a constant, under the condition that the maximum electrolyte packing concentration $n^0_s$ from other electrolytes is far less than the surfactant maximum packing concentration $\varGamma _{\infty }$ on the interface, we can derive the one-dimensional long-wave model similar as in Conroy et al. (Reference Conroy, Matar, Craster and Papageorgiou2011). The specific derivation process of the one-dimensional model is outlined in Appendix A.

By choosing the unperturbed radius of the thread $R_0$ as the length scale, $\sigma _0/\mu _1$ as the velocity scale, $R_0\mu _1/\sigma _0$ as the time scale, $\sigma _0/R_0$ as the pressure scale, $k_BT/e$ as the electric potential scale, $2\varGamma _{\infty }e N_A$ as scale for the charge density and $\varGamma _{\infty }$ as the scale for surfactant concentration, the final dimensionless one-dimensional equations are as follows:

(2.23) \begin{align}&\qquad\qquad\qquad q_t+q_z w+\frac {1}{2}qw_z=\frac {{\textit{Pe}}_s^{-1}}{S}\frac {\partial }{\partial z}(Sq_z+Sq\phi _{s,z}), \end{align}

(2.24) \begin{align}& \left [ \frac {S}{2}Re(w_t+{ww}_z)-3S_z w_z-\frac {3}{2}{\textit{S}w}_{\textit{zz}} \right ]\nonumber\\ &\qquad + \frac {S}{4}Q\left \{ \chi ^2\phi _s^2\left [ \frac {I_1^2(\chi S)}{I_0^2(\chi S)}-1 \right ] -\frac {\varepsilon _2}{\varepsilon _1}\left [ \frac {\phi _s-\phi _w}{S\ln (R_2/S)} \right ]^2 \right \}_z+2\beta q\phi _{s,z} \nonumber\\ &\quad =\left [ \sigma +q({\textit{Bo}}_{\kappa \infty }-{\textit{Bo}}_{\mu \infty })\left ( \frac {(\textit{S}w)_z}{S}+\zeta (-\frac {{\textit{S}w}_z}{2}-S_zw) \right ) \right ]_z \nonumber\\ &\qquad \quad \unicode{x2460}\quad \quad \quad \quad \quad \quad \quad \unicode{x2461} \nonumber\\ &\qquad -\frac {S}{2}\left \{ \zeta \left [ \sigma +q({\textit{Bo}}_{\kappa \infty }-{\textit{Bo}}_{\mu \infty })\left ( \frac {({\textit{S}w})_z}{S}+\zeta (-\frac {{\textit{S}w}_z}{2}-S_zw) \right ) \right ] \right \}_z \nonumber\\ &\qquad \quad \quad \quad \quad \unicode{x2462} \quad \quad \quad \quad \quad \unicode{x2463} \nonumber\\ &\qquad +\frac {5{\textit{Bo}}_{\mu \infty }}{2S}\boldsymbol{\cdot }\frac {\partial }{\partial z}[q{\textit{S}w}_z], \nonumber\\ &\quad \quad \quad \unicode{x2464} \end{align}

(2.25) \begin{align}&\qquad\qquad\qquad\qquad\qquad S_t+{w\textit{S}}_z+\frac {1}{2}{\textit{S}w}_z=0, \end{align}

where

(2.26) \begin{align}& \sigma =1+\beta \ln (1-2q)+\beta q\phi _s-\beta \frac {(\varepsilon _2/\varepsilon _1)\phi _w}{\chi ^2\xi S\ln (R_2/S)}(\phi _s-\varPhi ),\\[-12pt] \nonumber \end{align}

(2.27) \begin{align}&\quad\,\, \phi _s=I_0(\chi S)\frac {\chi ^2\xi qS\ln (R_2/S)+(\varepsilon _2/\varepsilon _1)\phi _w}{\chi I_1(\chi S)S\ln (R_2/S)+(\varepsilon _2/\varepsilon _1)I_0(\chi S)},\\[-12pt] \nonumber \end{align}

(2.28) \begin{align}&\qquad\qquad \varPhi =\frac {(\varepsilon _2/\varepsilon _1)\phi _w}{\varepsilon _2/\varepsilon _1+\chi \frac {I_1(\chi S)}{I_0(\chi S)}S\ln (R_2/S)},\\[-12pt] \nonumber \end{align}

(2.29) \begin{align}&\qquad\qquad\,\, \zeta =\frac {1}{S\sqrt {1+S_z^2}}-\frac {S_{\textit{zz}}}{\sqrt {(1+S_z^2)^3}}.\\[-12pt] \nonumber \end{align}

The one-dimensional model presented above comprises three independent variables: $q$ represents the surface charge density, which is related to the surfactant concentration, $q=( {1}/{2})\varGamma$ with initial values of $q_0$ and $\Gamma_0$ , respectively; $w$ represents the axial velocity, with an initial value of $0$ ; and $S$ indicates the position of the gas–liquid interface.

Equation (2.23) represents the transport equation for the surface charge, with the convection term on the left-hand side and the diffusion term on the right-hand side. Equation (2.24) is a one-dimensional momentum equation. The second and third terms on the left-hand side of the equation represent the electrostatic stress in the bulk stream and the electrostatic tension from the tangential stress equilibrium, respectively. On the right-hand side of this equation, term ① represents the Marangoni stress caused by the spatial inhomogeneity of surface tension and the electrical effect, term ③ is the capillary stress, and terms ②, ④ and ⑤ represent the surface viscous stresses associated with surface dilatational deformation and surface shear deformation.

The dimensionless numbers that appear in the above one-dimensional model are as follows: the Péclet number on the interface ${\textit{Pe}}_s=(\sigma _0 R_0)/(D_s \mu )$ , the Reynolds number ${\textit{Re}}=(\rho \sigma _0R_0)/(\mu ^2)$ , the electrostatic Weber number $Q=(\varepsilon _1 k_BT/e)/(\sigma _0 R_0)$ characterizing the ratio of electrostatic stress to capillary stress, the Marangoni number $\beta =R_g T\varGamma _{\infty }/\sigma _0$ quantifying the sensitivity of surface tension to surfactant concentration, the dilatational Boussinesq number ${\textit{Bo}}_{\kappa \infty }=(2\kappa _{\textit{sr}}\varGamma _{\infty })/(\mu _1 R_0 \varGamma _r)$ and the shear Boussinesq number ${\textit{Bo}}_{\mu \infty }=(2\mu _{\textit{sr}}\varGamma _{\infty })/(\mu _1 R_0 \varGamma _r)$ used to characterize the relative importance of surface viscosity versus bulk viscosity. Here, an additional parameter $\varTheta = \kappa _s/\mu _s={\textit{Bo}}_{\kappa \infty }/{\textit{Bo}}_{\mu \infty }$ is introduced to represent the dilatational-to-shear surface viscosity ratio. The remaining parameters in the one-dimensional model are as follows: $\chi =\sqrt {(2R_0^2 e^2 c_0 N_A)/(k_B T\varepsilon _1)}$ represents the ratio of the radius of the liquid thread to the length of Debye, $\xi =\varGamma _{\infty }/(c_0 R_0)$ is related to the concentration of surfactant and ions in the mainstream, and $\phi _s$ is the electric potential at the gas–liquid interface.

For an isothermal saturated interface, the Marangoni stress can be neglected, and the only relevant surface stresses are the Young–Laplace pressure and surface viscous stresses (Quéré & de Ryck Reference Quéré and de Ryck1998; Scheid et al. Reference Scheid, Delacotte, Dollet, Rio, Restagno, van Nierop, Cantat, Langevin and Stone2010; Scheid et al. Reference Scheid, Dorbolo, Arriaga and Rio2012; Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2020). Martinez–Calvo provided a detailed description of the applicable conditions, and as discussed by Quéré & de Ryck (Reference Quéré and de Ryck1998) in the context of axisymmetric flows, the validity of the assumption is ensured if two conditions are met: first, there must be sufficient surfactant available at the interface during the thinning process; second, the adsorption rate must exceed the characteristic velocity of the interface. This aligns with the diffusion-controlled model and the small Péclet number assumption used in this paper. Furthermore, when the surfactant concentration is low, the influence of the Marangoni effect can similarly be neglected (Ponce-Torres et al. Reference Ponce-Torres, Montanero, Herrada, Vega and Vega2017), as its impact on the surface tension coefficient is minimal. In addition, the surface viscosity coefficient is difficult to measure under significant Marangoni stress and adsorption dynamics, as all these effects are inherently entangled (Prosser & Franses Reference Prosser and Franses2001; Elfring, Leal & Squires Reference Elfring, Leal and Squires2016). Therefore, although the Marangoni effect was considered during the derivation process, it was neglected in the result calculations, which means $\beta =0$ . This decision is based on two reasons: first, under the aforementioned applicable conditions, neglecting the Marangoni effect is reasonable; second, more importantly, it serves to highlight the influence of surface viscosity, allowing for a more in-depth study of the mechanisms by which surface viscosity and electric effects influence the instability of liquid jets.

If the concentration of the surfactant is set to zero, that is to say the density of the surface charge is also set to zero, then the surface viscosity will also be absent. If the surface viscosity be set to zero, the one-dimensional model will degenerate into the model of Conroy et al. (Reference Conroy, Matar, Craster and Papageorgiou2011).

5.1. The evolution of thread radius over time

The effects of surface charge density $q_0$ and the two Boussinesq numbers characterizing surface viscosity on the evolution of the liquid thread over time are illustrated in figure 6. The curves represent the variation of the minimum radius of the liquid thread, $S_{\textit{min}}$ , over time.

Figure 6.

The variation of the radius at the thinnest point of the liquid thread over time. The effects of (a) the Boussinesq number and (b) the initial surface charge density. (c) The contributions of various stresses during the thinning process of the liquid thread with $q_0=0.075,{\textit{Bo}}_{\kappa \infty }={\textit{Bo}}_{\mu \infty }=10$ .

Figure 6(a) illustrates the effect of surface viscosity on the thinning dynamics of the liquid thread for $ q_0=0.075$ . Higher Boussinesq numbers significantly delay thinning, with a more pronounced effect at larger values. This result is consistent with the findings of Wee et al. (Reference Wee, Wagoner, Garg, Kamat and Basaran2021) and Li & He (Reference Li and He2023). Figure 6(b) examines the role of surface charge at $ {\textit{Bo}}_{\kappa \infty }={\textit{Bo}}_{\mu \infty }=10$ , showing that increased $ q_0$ reduces the thinning rate. This effect primarily arises from enhanced electrostatic repulsion, as the Marangoni effect is neglected here.

For small values of $ q_0$ and the Boussinesq numbers, the $ S_{\textit{min}}$ -time curve exhibits an inflection point when $ S_{\textit{min}} \lt 0.1$ , which disappears at larger values. Initially, the thinnest point of the liquid thread is located at $ L_{\lambda }/2$ . As nonlinear effects grow, this thinnest point shifts outward, slowing the thinning rate temporarily before transitioning to an exponential decay. In contrast, at higher $ q_0$ and Boussinesq numbers, the thinning remains centred at $ L_{\lambda }/2$ , making this effect less prominent. A detailed explanation follows in § 5.3.

Based on this foundation, the contributions of various stresses during the thinning process of the liquid thread have also been analysed. Considering the conditions $\beta =0$ and ${\textit{Bo}}_{\kappa \infty }={\textit{Bo}}_{\mu \infty }$ , (2.24) can be rewritten as follows:

(5.1) \begin{align} 0&=-\textit{Re}(w_t+{{ww}}_z)+(3w_{\textit{zz}}+\dfrac {6}{S}S_{z}w_{z})-\dfrac {1}{2}Q\left \{ \chi ^2\phi _s^2\left [ \frac {I_1^2(\chi S)}{I_0^2(\chi S)}-1 \right ]\right .\nonumber\\&\quad\left. -\frac {\varepsilon _2}{\varepsilon _1}\left [ \frac {\phi _s-\phi _w}{S\ln (R_2/S)} \right ]^2 \right \}_z -\left [ \zeta \sigma \right ]_z +\frac {5{\textit{Bo}}_{\mu \infty }}{S^2}\boldsymbol{\cdot }\frac {\partial }{\partial z}[q{\textit{S}w}_z]. \end{align}

The initial term on the right-hand side of (5.1) represents the inertial stress ( $f_I$ ), and indeed, the quantity represented by the term $(w_t + {{ww}}_z)$ characterizes the acceleration. The second term represents the viscous stress exerted by the liquid ( $f_l$ ), the third term represents the electrostatic stress ( $f_{{\kern-0.5pt}E}$ ), the fourth term represents the capillary stress ( $f_C$ ) and the fifth term represents the surface shear deformation viscous stress ( $f_S$ ). By substituting the calculated results into the expressions for each stress, one can determine the instantaneous magnitudes and directions of the stresses, with positive values indicating the positive direction along the z-axis and negative values indicating the opposite.

Figure 6(c) presents the stress balance during the initial and intermediate stages of liquid thread thinning. Capillary stress drives liquid outward, promoting thinning at the centre, while opposing stresses resist this process and inhibit breakup. Computational results confirm that capillary stress oppose other stresses along the $ z$ -axis, consistent with the analysis of Wee et al. (Reference Wee, Wagoner, Garg, Kamat and Basaran2021). Under the current parameters, surface shear and bulk viscous stress dominate this resistance, whereas electrostatic effects remain weak due to the low wall potential. A detailed analysis of cases with stronger electrostatic stress follows in subsequent sections.

5.2. The satellite droplets and main droplets

The formation of satellite droplets is a topic of great interest among researchers. This study examines the effects of the initial surface charge density $q_0$ and the Boussinesq numbers characterizing surface viscosity on the formation of satellite droplets.

The sizes of satellite and main droplets are characterized by their equivalent radius, $ R_s$ and $ R_m$ , respectively. The equivalent radius is defined as $ R_{s,m} = ( {3V_{s,m}}/{4\pi } )^{1/3}$ , where $ V_s$ and $ V_m$ denote the volumes of the satellite and main droplets. These volumes are obtained by integrating the interface shape at the moment of breakup.

Figure 7 illustrates the influence of $ q_0$ and the Boussinesq numbers on satellite droplet formation. Increasing $ q_0$ at fixed Boussinesq numbers (figure 7 a,b) or increasing Boussinesq numbers at fixed $ q_0$ (figure 7 a,c and 7 b,d) leads to smaller satellite droplets, consistent with the experimental findings of Ponce-Torres et al. (Reference Ponce-Torres, Montanero, Herrada, Vega and Vega2017), which demonstrated that surface viscosity reduces satellite droplet size. The details of a qualitative comparison with this experiment are presented in Appendix C. When both parameters are sufficiently large, satellite droplets nearly vanish (figure 7 d), suggesting a threshold in the $(q_0, {\textit{Bo}}_{\mu \infty })$ space beyond which satellite formation is suppressed.

Figure 7.

The effect of initial charge density and the Boussinesq numbers on the size of satellite and primary droplets at the breakup moment. (a–d) An increase in initial charge density and the Boussinesq numbers reduces the satellite droplet size while enlarging the primary droplet. (e) Phase diagram and boundary separating the presence and absence of satellite droplets. ( f) The equivalent radius of the primary droplet, $ R_m$ , exhibits a reliable correlation with the most unstable wavenumber, $ k_{\textit{max}}$ .

Figure 7(e) presents the phase diagram delineating the boundary for satellite droplet formation, defined by whether the thread radius at the midpoint exceeds $10^{-3}$ at the end of the simulation. Red inverted triangles denote cases without satellite droplets, blue circular markers indicate their presence and green squares represent the transitional boundary. This boundary is specific to the current parameter set and may shift with variations in other parameters, such as wall potential. The phase diagram reveals that increasing $ q_0$ and the Boussinesq numbers inhibits satellite droplet formation, primarily due to the reduced maximum growth rate. Although a small most unstable wavenumber typically favours satellite formation, the dominant effects of surface charge and viscosity counteract this trend.

In addition to satellite droplets, the size of the primary droplets is also a significant result in liquid stability studies since Lord Rayleigh’s seminal work (Rayleigh Reference Rayleigh1878).

Figure 7(a–d) illustrates that as $ q_0$ or the Boussinesq numbers increase, the main droplet radius $ R_m$ also increases. This is attributed to the decrease in the most unstable wavenumber, which leads to a longer dominant wavelength and, consequently, a larger primary droplet. However, the main droplet radius maintains a robust correlation with the wavenumber. As shown in figure 7( f), $k_{\textit{max}} {\times}R^3_m$ approaches a constant value of approximately $ 1.5\pi \approx 4.712$ , with the black dashed line representing the fitted curve. As $ k_{\textit{max}}$ increases, $k_{\textit{max}} {\times }R^3_m$ gradually declines, as lower surface charge and Boussinesq numbers promote satellite formation.

5.3. The mechanism by which surface viscosity suppresses satellite droplet formation

Section 5.2 shows that higher surface viscosity suppresses satellite droplet formation, warranting further investigation of the underlying mechanisms. We compare two cases: $ {\textit{Bo}}_{\kappa \infty } = {\textit{Bo}}_{\mu \infty } = 1$ (figure 8 a–e) and $ {\textit{Bo}}_{\kappa \infty } = {\textit{Bo}}_{\mu \infty } = 100$ (figure 8 g–n), both with $ q_0 = 0.075$ . Since electrostatic stress is negligible here (§ 5.1), the dominant stresses are capillary stress $ f_C$ and viscous stress $ f_V$ , with $ f_V = f_S + f_l$ , where $ f_S$ and $ f_l$ denote surface shear and bulk viscosity stress, respectively.

Figure 8.

Under a small Boussinesq number ${\textit{Bo}}_{\kappa \infty }={\textit{Bo}}_{\kappa \infty }=1$ , panel (a) shows the thinning process of the liquid thread, panel (b) illustrates the variation of velocity distribution over time, panel (c) presents the distribution of the liquid thread radius gradient $S_z$ at $t=54$ , panel (d) depicts the velocity gradient $w_z$ at $t=54$ and panel (e) displays the distribution of various stresses at $t=54$ . For a large Boussinesq number ${\textit{Bo}}_{\kappa \infty }={\textit{Bo}}_{\kappa \infty }=100$ , panel (g) depicts the thinning process of the liquid thread, panel (h) shows the velocity distribution over time, panel (l) presents the liquid thread radius gradient $S_z$ at $t=488$ , panel (m) illustrates the velocity gradient $w_z$ at $t=488$ and panel (n) displays the distribution of various stresses at $t=488$ . The initial charge density is $q_0=0.075$ for all cases.

For $ {\textit{Bo}}_{\kappa \infty } = {\textit{Bo}}_{\mu \infty } = 1$ , the breakup occurs at $ t_b = 57.6$ . Figure 8(a) depicts the evolution of the liquid thread for $ t \in [0,54]$ . At $ t = 54$ , the central section within a single wavelength forms an elongated liquid bridge. If this evolution continues, a satellite droplet will emerge at the midpoint, with the neck point forming between the main and satellite droplets. During thinning, capillary stress drives liquid outward, while viscosity resists this motion. As the bridge region ( $ z \in [-5,5]$ ) becomes thinner, viscous effects intensify, preventing complete liquid expulsion. Only the outer regions ( $ z \notin [-4,4]$ ) experience significant drainage, leading to liquid retention at the centre and an outward shift of the neck point. Figure 8(b) shows that the maximum velocity location (asterisk) moves outward as the thread thins, indicating that the site of maximum instantaneous liquid expulsion shifts outward, which drives the neck point’s migration. Near breakup ( $ t = 56.5$ ), velocity in the central region reverses, suggesting liquid accumulation at the midpoint to form a satellite droplet. From (5.1), the radius gradient $ S_z$ and velocity gradient $ w_z$ play a crucial role in balancing capillary and viscous stresses. Figure 8(c–e) show that stress peaks closely align with the peak of $ w_z$ , with a trend of stress reversal near the bridge centre.

For $ {\textit{Bo}}_{\kappa \infty } = {\textit{Bo}}_{\mu \infty } = 100$ , the breakup occurs at $ t_b = 488.6$ . Figure 8(g) shows the evolution of the liquid thread for $ t \in [0,488]$ , where the neck point remains at the centre throughout thinning, and no elongated liquid bridge forms. Without such a bridge, capillary stress continuously drives liquid outward, maintaining the neck point at the centre. From the velocity distribution in figure 8(h), the maximum velocity point shifts inward, accelerating liquid outflow near the centre and ensuring that the thinnest region remains at the midpoint. This behaviour contrasts sharply with the $ {\textit{Bo}}_{\mu \infty } = 1$ case, where the film shape near breakup exhibits steeper gradients, leading to bridge formation. The presence of a liquid bridge is a key determinant of satellite droplet formation. A higher surface viscosity suppresses surface gradients, as surfactant accumulation in thicker regions enhances local viscosity, leading to a smoother film profile that prevents bridge formation and inhibits satellite droplet generation. Figure 8(l) further confirms that the thickness gradient $ S_z$ is significantly lower for $ {\textit{Bo}}_{\mu \infty } = 100$ than for $ {\textit{Bo}}_{\mu \infty } = 1$ . The velocity gradient influences stress, which in turn affects the velocity distribution, driving the inward movement of the velocity peak.

It is worth noting that the above analysis is based on the influence of interfacial viscosity on the drainage rate and the morphology of the liquid thread, and it has already been stated at the beginning of this section that the effect of the electric field under the present conditions is weak. In § 5.4, we will further demonstrate that the decisive factor governing the distribution of surfactant concentration is the thread radius, while the electric field plays only a secondary role. This implies that the electric field cannot significantly modify the role of interfacial viscosity by altering the surfactant distribution. Therefore, the suppression of drainage and the inhibition of satellite droplet formation by interfacial viscosity can be regarded as intrinsic features, independent of the applied electric field. Consequently, it may be anticipated that the same interfacial-viscosity-driven mechanism would persist even in the presence of an axial electric field, which is not considered in the present work.

5.4. The mechanism of electrostatic stress and synergy with viscous stress

In the results presented in § 5.1, the electric stress was relatively small and did not significantly influence the breakup of the liquid thread. All aforementioned results were obtained with a wall potential of $\phi _w = 2$ . Therefore, to investigate how the results differ when the electric stress is comparable in magnitude to the viscous stress induced by surface shear deformation, it is worthwhile studying the effects of increasing the wall potential to enhance the electric stress.

In fact, Conroy et al. (Reference Conroy, Matar, Craster and Papageorgiou2011) indicated that as the wall potential increases, the electric stress rises, thereby enhancing the stability of the liquid thread. However, Conroy’s model did not account for the issue of surface viscosity. Here, by considering the surface viscosity introduced by surfactants, some different results have been observed.

The influencing factors of electric stress are first explored. The parameters are set as $ q_0 = 0.1$ and $ {\textit{Bo}}_{\kappa \infty } = {\textit{Bo}}_{\mu \infty } = 10$ . Figure 9(a) presents the effect of wall potential on the dispersion relation from linear stability results, showing that an increase in wall potential slightly reduces the maximum growth rate, though the variation is minimal, indicating a negligible influence. However, figure 9(b) reveals that a higher wall potential prolongs the final stage of liquid thread thinning, highlighting the significant role of electric stress in this phase. The study by Lakdawala et al. (Reference Lakdawala, Sharma and Thaokar2016) has also demonstrated that the influence of electric effects becomes significantly more pronounced in the final stage of liquid thread thinning compared with its early evolution.

Figure 9.

The influencing factors of electrostatic stress. (a) The effect of wall potential $\phi _w$ on the dispersion relation. (b) The influence of wall potential $\phi _w$ on the thinning process of the liquid thread. (c) The shape of the liquid thread near rupture under different wall potentials. (d) An enlarged view of the liquid thread shape approaching rupture. (e) The distributions of charge density, surface potential, surface potential gradient, the product $ \phi _s \boldsymbol{\cdot }\phi _{s,z}$ and the electrostatic stress. The parameter values are set as $q_0=0.1$ and ${\textit{Bo}}_{\kappa \infty }={\textit{Bo}}_{\kappa \infty }=10$ . We select the moment at $ t = 0.995 t_b$ , where $ t_b$ is the breakup time.

From (5.1), the electric stress $ f_{{\kern-0.5pt}E}$ depends on the product of the surface potential $ \phi _s$ and its gradient $ \phi _{s,z}$ . According to (2.27), $ \phi _s$ is negatively correlated with surface charge density $ q$ . In addition, from figure 9(d) and the charge distribution $q$ in figure 9(e), it is clear that the dominant factor controlling the charge concentration is the thread radius $S$ . Thus, the liquid thread radius influences the magnitude of electric stress, which in turn tends to draw liquid from the main droplet towards the origin. This effect becomes more pronounced as the wall potential increases. The results in figure 9(c–e) support these analyses. Figure 9(c–d) depict the liquid thread radius near rupture under different wall potentials, while figure 9(e) shows surface charge density, revealing a strong positive correlation with thread radius. The surface potential $ \phi _s$ is negatively correlated with $ q$ , while electric stress $ f_{{\kern-0.5pt}E}$ follows the trend of $ \phi _s \boldsymbol{\cdot }\phi _{s,z}$ .

To further analyse the electric potential and field distribution for $\phi _w = 4.5$ , the results are computed using (A35)–(A37) and presented in figure 10(a). The red curve marks the liquid thread interface. Although $R_2 = 5$ , only the region $r \lt 3$ is shown. The electric field outside the liquid thread is nearly perpendicular to the interface. Since the surfactant adsorbed on the surface carries a positive charge, it experiences an electric pressure that forces liquid outward at satellite droplets, leading to their reduction. This observation aligns with the findings of Conroy et al. (Reference Conroy, Matar, Craster and Papageorgiou2011). Inside the liquid thread, according to (A32) and the Debye–Hückel limit, the bulk charge is negative and related to $\phi _{10}$ . As depicted in figure 10(b), the internal electric field has a tangential component directed from the centre towards the ends, generating an electric stress $f_{{\kern-0.5pt}E}$ that drives fluid from both ends towards the centre, thereby suppressing liquid drainage. The key influencing factors of $f_{{\kern-0.5pt}E}$ have been detailed in figure 9.

Figure 10.

The mechanism of electric effects. (a) Distribution of electric potential and field at $\phi _w=4.5$ . (b) Local magnification of electric potential and field distribution, focusing on the internal electric field of the liquid thread. (c) Schematic model illustrating the twofold effects of electric stresses. Red arrows denote the electric field direction, blue arrows indicate internal electric stress.

Thus, the electric field exerts two primary effects: (i) externally, it applies a normal pressure on the surface-charged surfactant, reducing satellite droplet formation, and (ii) internally, it generates a converging electric stress on the negatively charged bulk, suppressing liquid drainage. The schematic in figure 10(c) illustrates these mechanisms. In addition to these electrical effects, capillary stress drives liquid drainage and satellite droplet formation, while surface viscosity counteracts the flow, inhibiting both drainage and satellite droplet formation. Additionally, it should be emphasized that the electric effect influences both the variation in thread thickness and the distribution of surfactants, while the surface charge density $ q$ is primarily related to the thread radius $ S$ , and surface viscosity is governed by surfactant concentration. Consequently, for the liquid region driven inward by the electric effect ( $ z \in [-8,-2] \cup [2,8]$ ), surface viscosity acts as a stabilizing factor.

Furthermore, ionic surfactants play a key role here by determining the type of surface charge. In conventional electrolyte solutions, the surface charge is governed by the wall potential: a positive wall potential leads to the accumulation of negative charges at the interface. In contrast, for the ionic surfactants considered here, whose adsorption characteristic time is much shorter than the diffusion characteristic time, the surface charge is determined by the type of ions associated with the surfactant (the study of Fauvel et al. (Reference Fauvel, Trybala, Tseluiko, Starov and Bandulasena2022) can also indirectly support this point). This represents a fundamental difference between ionic surfactants and ordinary electrolytes. A cationic surfactant is chosen specifically in this study to highlight this distinction. At the same time, for an anionic surfactant, similar results to those reported in this work can be obtained if the wall potential is correspondingly adjusted to a negative value.