2.1. Physical model and governing equations

Figure 1 shows a sketch of the swirling liquid jet, in which a viscous jet with radius $R_1$ is surrounded by ambient gas with radius $R_2$ under cylindrical coordinates of $(z, r,\theta )$. The liquid jet (with density $\rho _1$ and dynamical viscosity $\mu _1$) flows with velocity $\boldsymbol {U}^*_1(r)=U_1^*(r)\boldsymbol {e}_z+W_1^*(r)\boldsymbol {e}_\theta$ and the ambient gas (with density $\rho _2$ and dynamical viscosity $\mu _2$) evolves with velocity $\boldsymbol {U}_2^*(r)=U_2^*(r)\boldsymbol {e}_z+W_2^*(r)\boldsymbol {e}_\theta$, where $U_i^*(r)$ and $W_i^*(r)$ ($i=1, 2$) denote the axial and azimuthal velocity components, respectively. The surface tension of the liquid–gas interface is denoted by $\gamma$. To make the system dimensionless, the characteristic length, velocity, time and pressure are chosen as $R_1, W_1^*, R_1/W_1^*$ and $\rho _1W_1^{*2}$, respectively, where $W_1^*$ represents the azimuthal velocity at the jet interface. As a consequence, the dimensionless parameters involved in the theoretical model include the Reynolds number $Re=\rho _1W_1^*R_1/\mu _1$, the Weber number $We=\rho _1W_1^{*2}R_1/\gamma$, the density ratio $J=\rho _2/\rho _1$, the viscosity ratio $N=\mu _2/\mu _1$ and the radius ratio $a=R_2/R_1$. It is notable that the value of $a$ must be chosen much larger than unity to ensure a sufficient width of the gas environment.

Figure 1. Sketch of a viscous swirling liquid jet emerging into gas surroundings under cylindrical coordinates ($z, r, \theta$). The Rankine vortex velocity profile and the non-uniform velocity profile with the form of parabolic and error functions are utilized to approximate the azimuthal and the axial basic flows, respectively.

The flows of the swirling liquid jet and the ambient gas are governed by the dimensionless Navier–Stokes equations which consist of the continuity equation and the momentum equation, i.e.

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

(2.2)\begin{gather} \frac{\partial\boldsymbol{u}_i}{\partial t}+\boldsymbol{u}_i\boldsymbol{{\cdot}}\boldsymbol{\nabla}\boldsymbol{u}_i={-}\left(\frac{1}{J}\right)^{\delta_{2i}}\boldsymbol{\nabla} p_i+\left(\frac{N}{J}\right)^{\delta_{2i}}\frac{1}{Re}\nabla^2\boldsymbol{u}_i, \end{gather}

where $\boldsymbol {u}_i$ ($=(u_i,v_i,w_i)$) is the velocity component in cylindrical coordinates $(z, r, \theta )$, $p_i$ denotes the pressure, $\delta _{2i}$ is the Kronecker symbol and $i=1,2$ represent the liquid and the gas phases, respectively.

The corresponding boundary conditions are also given. At the symmetric axis $r=0$, the bounded condition and the single-valued condition should be satisfied as follows:

(2.3)\begin{equation} u_1<\infty,\quad v_1<\infty,\quad w_1<\infty,\quad \lim_{r\rightarrow 0}\frac{\partial \boldsymbol{u}_1}{\partial\theta}=0,\quad \lim_{r\rightarrow 0}\frac{\partial p_1}{\partial\theta}=0. \end{equation}

At the interface ($r=1+\eta$, where $\eta$ denotes a small perturbation of the interface position), the velocity continuity condition, the kinetic boundary condition and the dynamic boundary condition are given as

(2.4)\begin{gather} \boldsymbol{u}_1=\boldsymbol{u}_2, \end{gather}

(2.5)\begin{gather} v_1=\left(\frac{\partial}{\partial t}+\boldsymbol{u}_1\boldsymbol{{\cdot}}\boldsymbol{\nabla}\right)\eta, \end{gather}

(2.6)\begin{gather} (\boldsymbol{T}_2-\boldsymbol{T}_1)\boldsymbol{{\cdot}}\boldsymbol{n}-\frac{1}{We}(\boldsymbol{\nabla}\boldsymbol{{\cdot}}\boldsymbol{n})\boldsymbol{n}=0, \end{gather}

where $\boldsymbol {n}$ $(=(-({\partial \eta }/{\partial z}),1,-({\partial \eta }/{r\partial \theta })))$ represents the unit normal vector of the interface and $\boldsymbol {T}_1$ $(=-p_1\delta _{jl}\boldsymbol {e}_j \boldsymbol {e}_l+({1}/{Re})[\boldsymbol {\nabla }\boldsymbol {u}_1+(\boldsymbol {\nabla }\boldsymbol {u}_1)^T])$ and $\boldsymbol {T}_2$ $(=-p_2\delta _{jl}\boldsymbol {e}_j \boldsymbol {e}_l+({N}/{Re})[\boldsymbol {\nabla }\boldsymbol {u}_2+(\boldsymbol {\nabla }\boldsymbol {u}_2)^T])$ represent the hydrodynamic stress tensors of the liquid jet and the surrounding gas, respectively.

The boundary condition at $r=a$ is

(2.7)\begin{equation} u_2<\infty,\quad v_2<\infty,\quad w_2<\infty,\quad p_2<\infty. \end{equation}

2.2. Basic velocity profiles

To carry out the linear instability analysis, a basic flow state must be given first. In this work, the velocity profile at the azimuthal direction is established as the Rankine vortex, with the specific form

(2.8)\begin{gather} W_1^*(r)={\varOmega} r, \quad {\rm at} \ r \le R_1, \end{gather}

(2.9)\begin{gather} W_2^*(r)=\frac{{\varOmega} R_1}{r}, \quad {\rm at} \ r>R_1, \end{gather}

where $\varOmega$ is the rotary angular velocity of the liquid jet, as sketched in figure 1. The Rankine vortex describes the steady state of the swirling liquid jet and has been commonly utilized in previous theoretical studies (Moore & Saffman Reference Moore and Saffman1972; Loiseleux et al. Reference Loiseleux, Chomaz and Huerre1998; Gallaire & Chomaz Reference Gallaire and Chomaz2003a). The use of the Rankine vortex has also proved to be valid by experiments and numerical simulations (Billant et al. Reference Billant, Chomaz and Huerre1998; Sahu et al. Reference Sahu, Chetan, Mahato, Kar, Das and Lakkaraju2022). The dimensionless forms of the azimuthal velocities correspond to

(2.10)\begin{gather} W_1(r)=r, \quad {\rm at} \ r \le1, \end{gather}

(2.11)\begin{gather} W_2(r)=\frac{1}{r},\quad {\rm at} \ r>1. \end{gather}

As for the basic flow in the axial direction, we utilize the parabolic and error functions for the velocity profiles of the liquid jet and the gas surroundings, respectively. The parabolic function is known as the solution of pipe flow, which can be derived analytically to approximate the velocity profile of the liquid jet. The error function has been widely employed to construct the basic flow profiles of the jet or the mixing layers in previous studies (Yecko, Zaleski & Fullana Reference Yecko, Zaleski and Fullana2002; Otto, Rossi & Boeck Reference Otto, Rossi and Boeck2013; Matas, Delon & Cartellier Reference Matas, Delon and Cartellier2018). In our previous work (Li et al. Reference Li, Luo, Si and Xu2014), we successfully applied the error functions to model the surrounding air stream of an electrified liquid jet, in which a good agreement between theoretical predictions and experimental measurements was reached. The dimensional form of the axial velocity profile can be written as

(2.12)\begin{gather} U_1^*(r)=U_S+2\left[U_S-\frac{Q}{{\rm \pi} R_1^2}\right]\left[\frac{r^2}{R_1^2}-1\right], \quad {\rm at} \ r \le1, \end{gather}

(2.13)\begin{gather} U_2^*(r)=(U_\infty-U_S){\rm{erf}}\left[\frac{\sqrt{\rm \pi}}{2}\frac{4\mu_1\left(U_S-\dfrac{Q}{{\rm \pi} R_1^2}\right)}{\mu_2(U_\infty-U_S)}\left(\frac{r}{R_1}-1\right)\right]+U_S, \quad {\rm at}\ r>1, \end{gather}

where $U_S$ is the axial velocity at the surface and $U_\infty$ is the axial velocity at the outer boundary of the ambient gas (i.e. $r=a$), as sketched in figure 1, and $Q$ is the flow rate of the liquid jet. More details on the derivation of the velocity profile can be found in Appendix A. The dimensionless form of the axial velocity profile corresponds to

(2.14)\begin{gather} U_1(r)=V[V_S+2(V_S-1)(r^2-1)],\quad {\rm at} \ r \le1, \end{gather}

(2.15)\begin{gather} U_2(r)=V\left[(W-V_S) {\rm{erf}}\left[\frac{\sqrt{\rm \pi}}{2}\frac{4(V_S-1)}{N(W-V_S)}(r-1)\right]+V_S\right],\quad {\rm at} \ r>1, \end{gather}

where $V={\bar {U}_1}/{{\varOmega } R_1}$, $V_S={U_S}/{\bar {U}_1}$ and $W={U_\infty }/{\bar {U}_1}$. Here, $\bar {U}_1={Q}/{{\rm \pi} R_1^2}$ is the average axial velocity of the jet. To satisfy the real flow situations, the value of $V_S$ is directly related to that of $W$. Specifically, when the axial velocity of the jet is larger than that of the gas (i.e. $W<1$), the value of $V_S$ should be less than 1 to maintain the continuity of velocity at the interface. If the axial velocity of the jet is smaller than that of the gas (i.e. $W>1$), the value of $V_S$ should be larger than 1. When the jet velocity is equal to the gas velocity (i.e. $W=1$), the velocity profiles become uniform and thus $V_S=1$. It is notable that in figure 1, only the situation of $W<1$ is given for the sketch.

As the basic velocity profiles in the axial direction are non-uniform, axial shear stress exists both at the jet interface and in the fluid bulk. The strength of axial shear stress is in direct proportion to the derivatives of the axial velocities, with the form

(2.16)\begin{gather} \frac{{\rm d}U_1(r)}{{\rm d}r}=4V(V_S-1)r,\quad {\rm at} \ r \le1, \end{gather}

(2.17)\begin{gather} \frac{{\rm d}U_2(r)}{{\rm d}r}=4V(V_S-1)\exp\left({-\left[\frac{2\sqrt {\rm \pi}(V_S-1)}{N(W-V_S)}(r-1)\right]^2}\right),\quad {\rm at} \ r>1. \end{gather}

If $V_S=1$, the values of ${{\rm d}U_1(r)}/{{\rm d}r}$ and ${{\rm d}U_2(r)}/{{\rm d}r}$ are equal to zero invariably, and the liquid jet and the surrounding gas both flow with uniform velocity without axial shear stress. For the non-uniform velocity profiles with $V_S \ne 1$, the absolute values of the velocity derivatives become larger as $V_S$ gradually diverges from 1, indicating the enhancement of shear stress. A stronger axial shear stress can also be reached as $V$ increases. The value of $W$ does not affect the strength of axial shear stress but mainly modulates the thickness of the velocity boundary layer of gas surroundings. Specifically, as the value of $W$ gradually deviates from 1, the thickness of the gas velocity boundary layer becomes larger due to the increase of the velocity difference between the liquid jet and the gas surroundings.

2.3. Linear instability analysis

In theoretical analysis, the physical quantities ($\boldsymbol {u}_i, p_i$, where $i=1, 2$) are all divided into two parts: the basic state (i.e. $\boldsymbol {U}_i$ and $P_i$) and small disturbed quantities (denoted by $\tilde {\boldsymbol {u}}_i$ and $\widetilde {p_i}$). Utilizing these forms and ignoring the high-order nonlinear terms, the linear equations which govern the development of small perturbations are derived as

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

(2.19)\begin{gather} \frac{\partial\tilde{\boldsymbol{u}}_i}{\partial t}+\tilde{\boldsymbol{u}}_i\boldsymbol{{\cdot}}\boldsymbol{\nabla} \boldsymbol{U}_i+\boldsymbol{U}_i\boldsymbol{{\cdot}}\boldsymbol{\nabla}\tilde{\boldsymbol{u}}_i ={-}\left(\frac{1}{J}\right)^{\delta_{2i}}\boldsymbol{\nabla}\tilde{p}_i+\left(\frac{N}{J}\right)^{\delta_{2i}}\frac{1}{Re}\nabla^2\widetilde{\boldsymbol{u}_i}. \end{gather}

At the axisymmetric axis $r=0$,

(2.20)\begin{equation} \frac{\partial\tilde{v}_1}{\partial\theta}-\tilde{w}_1=0,\quad \tilde{v}_1+\frac{\partial\tilde{w}_1}{\partial\theta}=0,\quad \frac{\partial\tilde{u}_1}{\partial\theta}=0,\quad \frac{\partial\tilde{p}_1}{\partial\theta}=0. \end{equation}

At the jet interface $r=1+\eta$, the kinetic and dynamic boundary conditions can be written as

(2.21)\begin{gather} \tilde{u}_1+\frac{{\rm d}U_1}{{\rm d}r}\eta=\tilde{u}_2+\frac{{\rm d}U_2}{{\rm d}r}\eta, \end{gather}

(2.22)\begin{gather} \tilde{v}_1=\tilde{v}_2, \end{gather}

(2.23)\begin{gather} \tilde{w}_1+\frac{{\rm d}W_1}{{\rm d}r}\eta=\tilde{w}_2+\frac{{\rm d}W_2}{{\rm d}r}\eta, \end{gather}

(2.24)\begin{gather} \tilde{v}_1=\frac{\partial\eta}{\partial t}+U_1\frac{\partial\eta}{\partial z}+W_1\frac{\partial\eta}{r\partial\theta}, \end{gather}

(2.25)\begin{gather} \frac{1}{Re}\left(\frac{\partial\tilde{v}_1}{\partial z}+\frac{\partial\tilde{u}_1}{\partial r}+\frac{{\rm d}^2U_1}{{\rm d}r^2}\eta\right)=\frac{N}{Re}\left(\frac{\partial\tilde{v}_2}{\partial z}+\frac{\partial\tilde{u}_2}{\partial r}+\frac{{\rm d}^2U_2}{{\rm d}r^2}\eta\right), \end{gather}

(2.26)\begin{gather} \tilde{p}_1-\tilde{p}_2+\frac{{\rm d}P_1}{{\rm d}r}\eta-\frac{{\rm d}P_2}{{\rm d}r}\eta-\frac{2}{Re}\frac{\partial\tilde{v}_1}{\partial r}+\frac{2N}{Re}\frac{\partial\tilde{v}_2}{\partial r}+\frac{1}{We}\left(\eta+\frac{\partial^2\eta}{\partial \theta^2}+\frac{\partial^2\eta}{\partial z^2}\right)=0, \end{gather}

(2.27)$$\begin{gather} \frac{1}{Re}\left[\frac{\partial\tilde{v}_1}{\partial\theta}+\frac{\partial\tilde{w}_1}{\partial r}+\eta\frac{{\rm d}^2W_1}{{\rm d}r^2}-\left(\tilde{w}_1-\eta W_1+\frac{{\rm d}W_1}{{\rm d}r}\eta\right)\right]\nonumber\\ \qquad -\frac{N}{Re}\left[\frac{\partial\tilde{v}_2}{\partial\theta}+\frac{\partial\tilde{w}_2}{\partial r}+\eta\frac{{\rm d}^2W_2}{{\rm d}r^2}-\left(\tilde{w}_2-\eta W_2+\frac{{\rm d}W_2}{{\rm d}r}\eta\right)\right]=0. \end{gather}$$

At $r=a$, the boundary condition corresponds to

(2.28)\begin{equation} \tilde{u}_2=\tilde{v}_2=\tilde{w}_2=\tilde{p}_2=0. \end{equation}

The normal mode decomposition is employed to study the development of small perturbations, where the perturbations of velocity components, pressure and the displacement of jet interface $\eta$ are all expanded with the Fourier form, i.e.

(2.29)\begin{equation} (\tilde{u}_i,\tilde{v}_i,\tilde{w}_i,\tilde{p}_i,\eta)=(\hat{u}_i(r),\hat{v}_i(r),\hat{w}_i(r),\hat{p}_i(r),\hat{\eta}) \exp({{\rm{i}}(k z+m\theta-\omega t)}), \end{equation}

where $k$ is the dimensionless axial wavenumber, $m$ is the dimensionless azimuthal wavenumber and $\omega$ is the dimensionless perturbation frequency. In temporal instability analysis, $\omega$ is a complex number with the form $\omega =\omega _r+\rm {i}\omega _{i}$, $m$ is a non-negative integer and $k$ is a real number. The imaginary and real parts of $\omega$ (denoted by $\omega _i$ and $\omega _r$, respectively) stand for the temporal growth rate and the frequency of disturbance, respectively. If $\omega _i<0$ for all values of $k$, the jet is stable invariably. However, if $\omega _i>0$ for some values of $k$, the jet can be unstable. The jet is neutrally stable at $\omega _i=0$ as the initial disturbance neither grows nor decays with time. The variation of $m$ stands for different azimuthal modes. Mode $m=0$ represents the axisymmetric mode in which the jet evolves with axisymmetric disturbance, and $m\geq 1$ stands for the helical modes in which the jet develops under non-axisymmetric perturbations. Without loss of generality, this study mainly focuses on the growth of disturbance under azimuthal modes $m=0$–$5$.

Taking (2.29) into (2.18)–(2.28), the governing equations can be written as

(2.30)\begin{gather} \frac{{\rm d}\hat{v}_i}{{\rm d}r}+\frac{\hat{v}_i}{r}+\frac{{\rm i}m}{r}\hat{w}_i+{\rm i}k\hat{u}_i=0, \end{gather}

(2.31)\begin{gather}-\left(\frac{1}{J}\right)^{\delta_{2i}}{\rm i}k\hat{p}_i-\hat{v}_i\frac{{\rm d}U_i}{{\rm d}r}+\left(\frac{N}{J}\right)^{\delta_{2i}}\frac{1}{Re}\left\{\frac{{\rm d}^2\hat{u}_i}{{\rm d}r^2}+\frac{1}{r}\frac{{\rm d}\hat{u}_i}{{\rm d}r}-\left[k^2+\frac{m^2}{r^2}\right.\right.\nonumber\\ \left.\left.+\,{\rm i}\left(k U_i+\frac{W_im}{r}\right)\left(\frac{J}{N}\right)^{\delta_{2i}}Re\right]\hat{u}_i\right\}={-}{\rm i}\omega\hat{u}_i, \end{gather}

(2.32)\begin{gather}-\left(\frac{1}{J}\right)^{\delta_{2i}}\frac{{\rm d}\hat{p}_i}{{\rm d}r}+\left(-\frac{1}{Re}\frac{2{\rm i}m}{r^2}\left(\frac{N}{J}\right)^{\delta_{2i}}+\frac{2W_i}{r}\right)\hat{w}_i+\frac{1}{Re}\left(\frac{N}{J}\right)^{\delta_{2i}}\left\{\frac{{\rm d}^2\hat{v}_i}{{\rm d}r^2}+\frac{{\rm d}\hat{v}_i}{rdr}\right.\nonumber\\ \left.-\left[k^2+\frac{m^2+1}{r^2}+{\rm i}\left(k U_i+\frac{mW_i}{r}\right)Re\left(\frac{J}{N}\right)^{\delta_{2i}}\right]\hat{v}_i\right\}={-}{\rm i}\omega\hat{v}_i, \end{gather}

(2.33)\begin{gather} -\left(\frac{1}{J}\right)^{\delta_{2i}}\frac{{\rm i}m}{r}\hat{p}_i+\left(\frac{1}{Re}\frac{2{\rm i}m}{r^2}\left(\frac{N}{J}\right)^{\delta_{2i}}-\frac{W_i}{r}-\frac{{\rm d}W_i}{{\rm d}r}\right)\hat{v}_i+\frac{1}{Re}\left(\frac{N}{J}\right)^{\delta_{2i}}\left\{\frac{{\rm d}^2\hat{w}_i}{{\rm d}r^2}+\frac{1}{r}\frac{{\rm d}\hat{w}_i}{{\rm d}r}\right.\nonumber\\ \left.-\left[k^2+\frac{m^2+1}{r^2}+{\rm i}\left(k U_i+\frac{mW_i}{r}\right)Re\left(\frac{J}{N}\right)^{\delta_{2i}}\right]\hat{w}_i\right\}={-}{\rm i}\omega\hat{w}_i. \end{gather}

Accordingly, the boundary conditions become as follows:

(2.34)\begin{gather} \hat{w}_1=\hat{v}_1=\frac{{\rm d}\hat{u}_1}{{\rm d}r}=\frac{{\rm d}\hat{p}_1}{{\rm d}r}=0,\quad {\rm at} \ r=0,\quad \text{for} \ {m}=0, \end{gather}

(2.35)\begin{gather} \hat{u}_1=\hat{p}_1=0,\quad \hat{v}_1+i\hat{w}_1=0,\quad 2\frac{{\rm d}\hat{v}_1}{{\rm d}r}+{\rm i}\frac{{\rm d}\hat{w}_1}{{\rm d}r}=0,\quad {\rm at} \ r=0, \quad {\rm for} \ {m}=1, \end{gather}

(2.36)\begin{gather} \hat{u}_1=\hat{v}_1=\hat{w}_1=\hat{p}_1=0, \quad {\rm at} \ r=0,\quad {\rm for} \ {m}\geq 2. \end{gather}

At the jet interface $r=1+\eta$,

(2.37)\begin{gather} \hat{u}_1+\frac{{\rm d}U_1}{{\rm d}r}\hat{\eta}=\hat{u}_2+\frac{{\rm d}U_2}{{\rm d}r}\hat{\eta}, \end{gather}

(2.38)\begin{gather} -{\rm i}\omega\hat{\eta}=\hat{v}_1-{\rm i}\left(k U_1+\frac{W_1 m}{r}\right)\hat{\eta}, \end{gather}

(2.39)\begin{gather} {\rm i}k\hat{v}_1+\frac{{\rm d}\hat{u}_1}{{\rm d}r}+\frac{{\rm d}^2U_1}{{\rm d}r^2}\hat{\eta}=N\left({\rm i}k\hat{v}_2+\frac{{\rm d}\hat{u}_2}{{\rm d}r}+\frac{{\rm d}^2U_2}{{\rm d}r^2}\hat{\eta}\right), \end{gather}

(2.40)\begin{gather} \hat{p}_1-\hat{p}_2+\frac{{\rm d}P_1}{{\rm d}r}\hat{\eta}-\frac{{\rm d}P_2}{{\rm d}r}\hat{\eta}-\frac{2}{Re}\frac{{\rm d}\hat{v}_1}{{\rm d}r}+\frac{2N}{Re}\frac{{\rm d}\hat{v}_2}{{\rm d}r}+\frac{1}{We}(1-m^2-k^2)\hat{\eta}=0, \end{gather}

(2.41)\begin{gather} {\rm i}m\hat{v}_1+\frac{{\rm d}\hat{w}_1}{{\rm d}r}-\hat{w}_1+\left(\frac{{\rm d}^2W_1}{{\rm d}r^2}+W_1-\frac{{\rm d}W_1}{{\rm d}r}\right)\hat{\eta}\nonumber\\ -\,N\left[{\rm i}m\hat{v}_2+\frac{{\rm d}\hat{w}_2}{{\rm d}r}-\hat{w}_2+\left(\frac{{\rm d}^2 W_2}{{\rm d}r^2}+W_2-\frac{{\rm d}W_2}{{\rm d}r}\right)\hat{\eta}\right]=0. \end{gather}

At the outside boundary $r=a$,

(2.42)\begin{equation} \hat{u}_2=\hat{v}_2=\hat{w}_2=\hat{p}_2=0. \end{equation}

2.4. Numerical solutions and validations

The non-uniform characteristics of the velocity profiles bring about a difficulty in deducing an analytical dispersion relation of perturbations. Therefore, we utilize the Chebyshev collocation method (Weideman & Reddy Reference Weideman and Reddy2000; Schmid & Henningson Reference Schmid and Henningson2001; Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021) to solve the perturbation growth numerically. In the process of solution, the liquid region $r\in [0,1]$ is mapped into the computational space of $y\in [-1,1]$ through the linear transformation

(2.43)\begin{equation} r=\frac{1+y}{2} \end{equation}

and the gas region $r\in [1,a]$ is mapped into the computational space $y\in [-1,1]$ by means of linear transformation

(2.44)\begin{equation} r=\frac{y(1-a)+(1+a)}{2}. \end{equation}

In this way, the physical quantities $\hat {u}_i(r)$, $\hat {v}_i(r)$, $\hat {w}_i(r)$, $\hat {p}_i(r)$ and their derivatives are expanded in Chebyshev polynomials and inserted into governing equations as well as boundary conditions to obtain discrete equations. The Gauss–Lobatto collocation points $y_j=\cos (j{\rm \pi} /N)$, $j=0,1,\ldots,N$, are employed to discretize the computation domain. Therefore, the discrete equations and the corresponding boundary conditions form a generalized eigenvalue problem with the form

(2.45)\begin{equation} {\rm i}\omega[\boldsymbol{B}]\boldsymbol{X}=[\boldsymbol{A}]\boldsymbol{X}, \end{equation}

where $\boldsymbol {X}$ denotes the characteristic vector. As the numbers of collocations for the liquid and the gas domains are $N_1+1$ and $N_2+1$, respectively, the sizes of two coefficient matrices $\boldsymbol {A}$, $\boldsymbol {B}$ are equal to $(4N_1+4N_2+9)\times (4N_1+4N_2+9)$. A MATLAB code is utilized to solve the generalized eigenvalue problem with QZ algorithm. In our study, we choose the value of $a=10$ to ensure the infinite characteristics of the outside boundary. The numbers of the collocation points are selected to be $N_1=30$ and $N_2=60$ to satisfy the calculation accuracy. The convergence studies on the selections of values of $a$, $N_1$ and $N_2$ and the numerical validations of our results to previous studies are given in Appendix B.

2.5. Energy budget

To analyse the mechanism of the jet instability, the energy budget equation is also built, which is able to trace the source of disturbance kinetic energy (Lin & Chen Reference Lin and Chen1998; Lin Reference Lin2003; Ding et al. Reference Ding, Mu, Si and Jian2022). Considering a control volume of the liquid jet over one perturbation wavelength $\lambda$, we form a dot product of (2.19) with dimensionless perturbation velocity $\tilde {\boldsymbol {u}}_1$ and use the continuity equation of (2.18) to simplify the pressure term and the Gauss theorem to reduce some of the volume integrals to surface integrals. Through integrating over the control volume and averaging over one wavelength $\lambda$ and one perturbation period $T=2{\rm \pi} /\omega _r$, the energy equation can be derived as

(2.46)$$\begin{gather} \frac{1}{T\lambda}\int_0^T\int_V\left(\frac{\partial}{\partial t}+\boldsymbol{U}_1\boldsymbol{{\cdot}}\boldsymbol{\nabla}\right)e\,{\rm d}V\,{\rm d}t={-}\frac{1}{T\lambda}\int_0^T\int_V\tilde{\boldsymbol{u}}_1\boldsymbol{{\cdot}}(\tilde{\boldsymbol{u}}_1\boldsymbol{{\cdot}}\boldsymbol{\nabla}\boldsymbol{U}_1)\,{\rm d}V\,{\rm d}t \nonumber\\ -\frac{1}{T\lambda}\int_0^T\int_A\tilde{p}_1\tilde{\boldsymbol{u}}_1\boldsymbol{{\cdot}}\boldsymbol{n}\,{\rm d}A\,{\rm d}t+\frac{1}{T\lambda Re}\int_0^T\int_A\tilde{\boldsymbol{u}}_1\boldsymbol{{\cdot}}\boldsymbol{\tau}_1\,{\rm d}A\,{\rm d}t-\frac{1}{2ReT\lambda}\int_0^T\int_V\boldsymbol{\tau}_1:\boldsymbol{\tau}_1\,{\rm d}V\,{\rm d}t, \end{gather}$$

where $e=\tfrac {1}{2}\tilde {\boldsymbol {u}}_1\boldsymbol {{\cdot }}\tilde {\boldsymbol {u}}_1$ denotes the disturbance kinetic energy and $V$ and $A$ stand for the control volume and surface area, respectively. Considering the tangential and normal dynamic boundary conditions at the jet surface, the energy budget can be finally obtained as

(2.47)$$\begin{align} \overline{KE}&=\overline{REY}+\overline{PRL}+\overline{SHL}+\overline{NVL}+\overline{PRG}+\overline{SUT}+\overline{NVG}\nonumber\\ &\quad +\,\overline{SHG}+\overline{AHG}+\overline{SHB}+\overline{AHB}+\overline{LSS}+\overline{OGS}+\overline{DIS}. \end{align}$$

The detailed expressions of each term in (2.47) are given in Appendix C. The left-hand-side term $\overline {KE}$ represents the time rate of change of the disturbance kinetic energy. For the right-hand-side terms, $\overline {REY}$ represents the energy transfer between the disturbance and the basic flow through the Reynolds stress; $\overline {PRL}$, $\overline {SHL}$ and $\overline {NVL}$ demonstrate the rate of work done by the pressure, the tangential and normal components of the viscous stress at the top and the bottom ends of the control volume, respectively. Since the temporal analysis assumes the jet is infinite and periodic in the axial direction, the terms worked at the top and the bottom ends of the control volume equal zero invariably. Therefore, $\overline {PRL}$, $\overline {SHL}$ and $\overline {NVL}$ are neglected in energy budget analysis. Term $\overline {SUT}$ is the rate of work done by the surface tension. Term $\overline {PRG}$ is the rate of work done by the gas pressure perturbation on the swirling jet. Term $\overline {NVG}$ represents the rate of work done by the normal viscous stress exerted by the perturbed gas at the interface. Terms $\overline {SHG}$ and $\overline {AHG}$ are the rate of work done by axial shear stress and azimuthal shear stress exerted by the perturbed gas at the interface, respectively. Terms $\overline {SHB}$ and $\overline {AHB}$ represent the rate of work done by the axial shear stress and azimuthal shear stress associated with the basic flow distortion caused by the interface displacement, respectively. Terms $\overline {LSS}$ and $\overline {OGS}$ stand for the work done by centrifugal force caused by the swirling of liquid and gas phases, respectively. Term $\overline {DIS}$ is the rate of mechanical energy dissipation through viscosity.

To eliminate the non-uniqueness of the energy terms, we utilize the disturbance kinetic energy $EK$ to normalize each term in (2.47), i.e.

(2.48)\begin{equation} KE=\overline{KE}/EK,\quad REY=\overline{REY}/EK,\quad PRL=\overline{PRL}/EK,\quad SHL=\overline{SHL}/EK,\quad{\rm etc.}, \end{equation}

with the specific form

(2.49)\begin{equation} EK=\int_0^T\int_0^{\lambda}\int_0^2\int_0^1\frac{\rm \pi}{2}(\tilde{u}_1^2+\tilde{v}_1^2+\tilde{w}_1^2)r\,{\rm d}r\,{\rm d}\theta \,{\rm d}z\,{\rm d}t. \end{equation}

This strategy has been widely used in previous studies (Ye, Yang & Fu Reference Ye, Yang and Fu2016; Ding et al. Reference Ding, Mu, Si and Jian2022). It is notable that a positive value of a certain term indicates that the corresponding force would promote the jet instability, while a negative value suggests that the corresponding force would suppress the instability.

2.6. Reference state

To investigate the effect of circumferential rotation and axial flow on jet stability, a reference state should be given first. To ensure the rationality of theoretical analysis, the reference state must correspond to the real flow situation. Following the previous experimental study of Kubitschek & Weidman (Reference Kubitschek and Weidman2007b), we choose silicone oil and air for the liquid jet and the surrounding gas, respectively. The densities and dynamic viscosities of fluids correspond to $\rho _1=971\, {\rm kg}\,{\rm m}^{-3}$, $\rho _2=1.29\, {\rm kg}\,{\rm m}^{-3}$, and $\mu _1=0.1223\,{\rm Pa}\,{\rm s}$, $\mu _2=2\times 10^{-5}\,{\rm Pa}\,{\rm s}$, respectively. The interfacial tension is $\gamma =26\,{\rm mN}\,{\rm m}^{-1}$ and the radius of the liquid jet is $R_1=3.5\,{\rm mm}$. For the reference state, the average axial velocity of the jet ($\bar {U}_1$) is chosen to be $2.6\,{\rm m}\,{\rm s}^{-1}$ and the annular velocity ($\varOmega$) is equal to $211\,{\rm rad}\,{\rm s}^{-1}$. The jet is supposed to move in static gas surroundings (i.e. $U_\infty =0$). We also assume the axial velocity $V_S$ on the interface is slightly slower than the average axial velocity of the jet. The specific values of the dimensionless parameters under the reference state are

(2.50)\begin{align} &Re=20.5, \quad We=71.3, \quad V=3.5, \quad V_S=0.995,\quad W=0,\quad J=0.0013,\nonumber\\ &\quad N=1.67\times 10^{{-}4}. \end{align}

It is notable that as the characteristic velocity scale is chosen as $\varOmega R_1$, the change of the rotating angular velocity would affect the values of $Re$ and $We$ simutaneously. Moreover, the value of $V$ ($=\overline {U_1}/\varOmega R_1$) also changes with the variation of $\varOmega$ in order to maintain a constant axial velocity $\overline {U_1}$. For the convenience of analysis, we also define the Ohnesorge number $Oh=\sqrt {We}/Re=\mu _1/\sqrt {\rho _1\gamma R_1}$ and the Reynolds number based on the axial velocity in the $z$ direction, i.e. $Re_z=Re V=\rho _1\bar {U}_1 R_1/\mu _1$. For the reference state, the value of $Oh$ is equal to 0.41, and the variations of angular velocity and axial velocity can be represented through singly changing the value of $Re$ and $Re_z$, respectively.

4.1. Effect of axial velocity profile

We firstly consider the effect of axial velocity profile on jet instability through changing the value of $V_S$. For the reference state, the ambient gas remains static with $W=0$, indicating that only the non-uniform velocity profiles with $V_S<1$ would conform to real flow situations. In our study, the value of $V_S$ is varied within the range from 0.9 to 1 to modulate the degree of shear stress, where $V_S=1$ represents the case of a discontinuous velocity profile (i.e. $U_1=V$ and $U_2=0$) without axial shear stress on the jet interface.

As we mainly focus on the mode transitions between different azimuthal wavenumbers $m$, we no longer give the growth rate curves and only display the maximum growth rates $\omega _{imax}$ accompanied by the most unstable wavenumbers $k_{max}$ under different values of $V_S$, as shown in figures 6(a) and 6(b), respectively. It can be seen from figure 6(a) that the values of $\omega _{imax}$ increase for all azimuthal modes as $V_S$ decreases from 1, suggesting that the instabilities of the perturbations will be strengthened as the axial shear stress is enhanced. Comparing the values of $\omega _{imax}$ between different azimuthal modes, the predominant mode can be determined, as sketched by the red numbers in figure 6(a). It is observed that the mode with smaller $m$ gradually becomes the predominant mode as $V_S$ decreases, suggesting that the increase of axial shear stress can promote the mode transition to smaller azimuthal wavenumber. It should be emphasized that the helical mode $m=1$ appears to be the predominant mode at $V_S=0.98$, which is not observed in § 3. The emergence of $m=1$ mode is mainly attributed to the non-negligible axial shear effect. Figure 6(b) shows that the decrease of $V_S$ also leads to an increase of $k_{max}$ for all azimuthal modes. The most unstable wavenumbers of the predominant mode under different $V_S$ are highlighted with filled symbols and connected with the black thick dashed line. It can be clearly seen that as the value of $V_S$ decreases, $k_{max}$ of the predominant mode increases monotonically, indicating that the enhancement of the axial shear stress could result in smaller size of droplets. Moreover, the values of $k_{max}$ remain significantly smaller than unity under relatively large $V_S$ (e.g. $V_S>0.98$) but increase to much larger than unity as $V_S$ decreases (e.g. $V_S \le 0.98$) for the predominant mode, showing that the jet breakup regime gradually shifts from the Rayleigh mode to the Taylor mode with an increase of shear stress, which is in accordance with previous work on a non-swirling liquid jet (Lin & Chen Reference Lin and Chen1998).

Figure 6. (a) Maximum perturbation growth rate $\omega _{imax}$ for different azimuthal modes as $V_S$ varies. (b) Most unstable wavenumber $k_{max}$ for different azimuthal modes as $V_S$ varies, where the most unstable wavenumbers corresponding to the predominant mode are depicted by the filled symbols and connected by the black thick dashed line. In both panels, the predominant modes are indicated by red numbers. The parameters are unified at $Re=20.5, Oh=0.41, V=3.5, W=0, J=0.0013$ and $N=1.67\times 10^{-4}$.

Figure 7 depicts the energy budget of the most unstable axial wavenumber $k_{max}$ as $V_S$ gradually changes. Similar to the analysis in figure 5, only the energy budgets of the predominant modes (indicated by the red numbers) are given in figure 7. Clearly, the terms $LSS$, $PRG$ and $SHG$ are positive, indicating that the centrifugal force of jet swirling, the gas pressure perturbation and the axial shear stress exerted by the gas surroundings at the interface can promote the growth of disturbances. As the terms $SUT$ and $DIS$ are negative, the surface tension and the viscous dissipation would restrain the perturbation growth. For the discontinuous velocity profile ($V_S=1$), the instability of the swirling liquid jet is mainly induced by the centrifugal force as the term $LSS$ has the maximum positive value. As $V_S$ gradually decreases, the terms $SHG$ and $PRG$ become larger simultaneously, which indicates that both the axial shear stress and the perturbation of gas pressure at the interface play a more significant role in jet instability, and the absolute value of $SHG$ remains larger than that of $PRG$. Since there is little difference in the energy budget results between the cases of $V_S=0.995$ and 1, the physical mechanism of jet breakup could hardly be affected under relatively weak axial shear stress. However, as $V_S$ keeps decreasing, the terms $SHG$ and $PRG$ gradually exceed $LSS$, indicating that the axial KHI is a main contribution for the growth of perturbations when the shear stress is relatively strong. In these situations, the jet instability presents a lower azimuthal wavenumber and finally converts to the axisymmetric mode with $m=0$. For the negative terms which suppress the jet instability, as $DIS$ has the smallest value, the viscous dissipation plays the leading effect.

Figure 7. Energy budget of the most unstable perturbation wave at the predominant mode as $V_S$ changes under $Re=20.5, Oh=0.41, V=3.5, W=0, J=0.0013$ and $N=1.67\times 10^{-4}$, where the predominant modes are indicated by red numbers.

4.2. Effect of axial velocity of liquid jet

The value of $V$ reflects the axial velocity of the swirling liquid jet. A larger $V$ under constant $V_S$ can result in faster axial velocity of the jet and thus stronger axial shear stress on the jet interface, as shown by (2.16) and (2.17), respectively. Figures 8(a) and 8(b) show the maximum perturbation growth rate $\omega _{imax}$ and the most unstable wavenumber $k_{max}$ for azimuthal modes for $m=0$–$5$ as the value of $V$ gradually varies, respectively. Clearly, the values of $\omega _{imax}$ for each mode increase gradually as $V$ increases, as shown in figure 8(a), indicating that the increase of axial velocity can promote the jet instability. It is also observed that the azimuthal mode with a smaller $m$ presents a faster increasing tendency of $\omega _{imax}$, which demonstrates that the mode with a smaller $m$ would be affected more strongly by the axial shear stress. Comparing the maximum value of $\omega _{imax}$ between different azimuthal modes, the increase of $V$ (i.e. the enhancement of shear stress) can lead to the predominant mode transition to smaller $m$, as indicated by the red numbers. Figure 8(b) shows an overall increasing tendency of $k_{max}$ for each azimuthal mode as $V$ increases. Moreover, the most unstable wavenumbers $k_{max}$ corresponding to the predominant mode are highlighted with filled symbols and connected with the black thick dashed line. It can be clearly seen that an increase of $V$ can lead to an increase of $k_{max}$ of the predominant mode, implying that the smaller size of droplets will be generated as the axial jet velocity increases. Similar to the situation of $V_S$ decreasing (see figure 6b), an enhancement of shear stress by increasing $V$ can also lead to the transition of jet breakup regimes from the Rayleigh mode to the Taylor mode.

Figure 8. (a) Maximum perturbation growth rate $\omega _{imax}$ for different azimuthal modes as $V$ varies. (b) Most unstable wavenumber $k_{max}$ for different azimuthal modes as $V$ varies, where the most unstable wavenumbers corresponding to the predominant mode at different $V$ are depicted by the filled symbols and connected by the black thick dashed line. In both panels, the predominant modes are indicated by red numbers. The parameters are unified at $Re=20.5, Oh=0.41, V_S=0.995, W=0, J=0.0013$ and $N=1.67\times 10^{-4}$.

The energy budget results of the predominant mode (see the red numbers) as $V$ gradually changes are shown in figure 9. Similarly, we only focus on the energy budget at the most unstable wavenumber $k_{max}$. The terms $LSS$, $PRG$ and $SHG$ are positive, indicating that the centrifugal force, the gas pressure perturbation and the axial shear stress exerted by the gas surroundings at the jet interface would promote the growth of disturbances. As the terms $SUT$ and $DIS$ are negative, the surface tension and the viscous dissipation would restrain the perturbation growth. Different from the situation of $V_S$ varying (see figure 7), the values of $PRG$ are larger than those of $SHG$ invariably, indicating that the gas pressure perturbation has a more significant effect than the shear stress on jet instability. With an increase of $V$, the values of $PRG$ and $SHG$ increase rapidly while the values of $LSS$ present a firstly increasing and then decreasing tendency. Specifically, $LSS$ remains the predominant energy term as $V$ is less than a critical value around 10, suggesting that the centrifugal force plays the primary role in jet instability at relatively low axial velocity. Therefore, the breakup of the liquid jet presents the helical modes with $m \ge 1$. As $V$ exceeds the critical value, $PRG$ becomes the predominant energy term, which indicates that the gas pressure perturbation at the interface mainly contributes to jet breakup under large axial velocity. The enhancement of KHI leads to the final transition of jet breakup to the axisymmetric mode with $m=0$. The value of $DIS$ remains negative with the smallest value among all energy terms, which indicates that the viscous dissipation plays the primary role in stabilizing the jet.

Figure 9. Energy budget of the most unstable perturbation wave at the predominant mode as $V$ gradually changes under $Re=20.5, Oh=0.41, V_S=0.995, W=0, J=0.0013$ and $N=1.67\times 10^{-4}$, where the predominant modes are indicated by red numbers.

4.3. Effect of axial velocity of gas surroundings

The parameter $W$ represents the ratio of average axial velocity between the ambient gas and the liquid jet. As $W$ gradually diverges from 1, the velocity boundary layer of the gas widens but the degree of shear stress at the interface maintains unchanged, as shown by (2.16) and (2.17). It has been indicated in § 2.2 that if $W>1$, the value of $V_S$ should also be larger than 1. Therefore, we choose $V_S=1.005$ for the cases of $W>1$ to maintain the same shear stress as the reference state (i.e. $W=0, V_S=0.995$). Also, we set $V_S=1$ as $W=1$ to realize equal velocity between the liquid and the gas, and the shear stress is absent.

The maximum perturbation growth rates $\omega _{imax}$ under different azimuthal modes as $W$ increases from 0 to 6 are presented in figure 10(a). As the axial shear effect is totally ignored at $W=1$, the value of $\omega _{imax}$ reaches the smallest. As $W$ diverges from 1, the values of $\omega _{imax}$ increase generally for each mode. Specifically, the maximum growth rates under $W=0$ and $W=2$ almost coincide with each other. It is also found that the modes with smaller azimuthal wavenumbers are more sensitive to the increase of $W$ as a faster increasing tendency of $\omega _{imax}$ can be observed. The red numbers in figure 10(a) demonstrate the predominant modes. When the value of $W$ gradually diverges from 1, the mode with smaller azimuthal wavenumber $m$ becomes the predominant mode, indicating that a larger velocity difference between the ambient gas and the liquid jet can promote jet breakup with smaller $m$. The most unstable wavenumbers $k_{max}$ for different azimuthal modes as $W$ varies are displayed in figure 10(b). It is found that $k_{max}$ of all modes increase generally as the value of $W$ gradually diverges from 1. The most unstable wavenumbers $k_{max}$ corresponding to the predominant mode are highlighted with filled symbols and connected with the black thick dashed line. It can be observed that $k_{max}$ of the predominant mode increases when the value of $W$ diverges from 1, indicating that a greater gas boundary layer thickness can lead to smaller droplet size and also promote the transition of jet instability regimes from Rayleigh mode to Taylor mode.

Figure 10. (a) Maximum perturbation growth rate $\omega _{imax}$ for different azimuthal modes as $W$ varies. (b) Most unstable wavenumber $k_{max}$ for different azimuthal modes as $W$ varies, where the most unstable wavenumbers corresponding to the predominant mode are depicted by the filled symbols and connected by the black thick dashed line. In both panels, the predominant modes are indicated by red numbers. The parameters are unified at $Re=20.5, Oh=0.41, V_S=0.995, V=3.5, J=0.0013$ and $N=1.67\times 10^{-4}$.

Figure 11 presents the energy budget of the predominant mode as $W$ gradually changes, in which the predominant modes are indicated by red numbers. Only the results corresponding to the most unstable wavenumber $k_{max}$ are given. Similar to figures 7 and 9, the values of $LSS$, $PRG$ and $SHG$ remain positive, while those of $SUT$ and $DIS$ are negative. It is notable that the term $PRG$ is markedly larger than $SHG$, especially under large values of $W$, indicating that the change of $W$ mainly affects the jet instability through gas pressure perturbation instead of shear stress. The physical mechanism lies in that the change of $W$ mainly modulates the thickness of the gas velocity boundary layer but does not change the degree of shear stress at the interface. Under relatively low gas velocity (e.g. $0\leq W\leq 3$), the change of $W$ has little impact on the values of all the energy terms, and the predominant mode remains at $m=3$. Once $W \ge 3$, the increase of $W$ leads to obvious changes of the energy terms $LSS$, $PRG$, $SUT$ and $DIS$. Specifically, the value of $LSS$ presents an overall decreasing tendency, while the value of $PRG$ increases invariably. As the predominant energy term finally changes from $LSS$ to $PRG$, the enhancement of gas pressure perturbation promotes the transition from helical mode to axisymmetric mode (i.e. $m=3$ to 0). For the negative energy terms, $DIS$ maintains the smallest value, and thus the viscous dissipation has the primary stabilizing effect.

Figure 11. Energy budget of the most unstable perturbation wave at the predominant mode as $W$ gradually changes under $Re=20.5, Oh=0.41, V_S=0.995, V=3.5, J=0.0013$ and $N=1.67\times 10^{-4}$, where the predominant modes are indicated by red numbers.