2. Methodology

In the following, we describe a conceptual framework for the analysis of linearly amplifying distortions in a liquid jet injected into a quiescent gas medium. The formulation is based on spatial stability analysis and accounts for the effects of viscosity as well as surface tension under the assumption of an axisymmetric laminar base state.

2.1. Base flow

The analyses employed in this work consider a linearization of the governing equations around a two-fluid base flow describing a liquid jet injected into quiescent gas. In general, the velocity distribution of axisymmetric liquid jets is governed by the distance to the nozzle exit, with the flow in the vicinity of the nozzle being influenced by the profile attained at the orifice. In an idealized configuration of a pipe with circular cross-section and constant radius connecting the nozzle to a reservoir, the product of the pipe radius $a$ and the Reynolds number introduces a characteristic length, $L_c=a\,Re$. If the ratio of the characteristic length to the length of the pipe is large, the profile at the nozzle exit is effectively uniform, leading to an intense shear layer at the gas–liquid interface past the nozzle exit. On the other hand, for long pipes with a comparatively small characteristic length, the flow at the nozzle exit effectively describes a Hagen–Poiseuille profile (Duda & Vrentas Reference Duda and Vrentas1967). Past the nozzle, the liquid begins to relax towards a uniform velocity distribution (González-Mendizabal, Olivera-Fuentes & Guzmán Reference González-Mendizabal, Olivera-Fuentes and Guzmán1987; Shkadov & Sisoev Reference Shkadov and Sisoev1996), with the region of influence of the upstream flow again scaling with the product of nozzle diameter and Reynolds number. For the case of a planar jet, Söderberg & Alfredsson (Reference Söderberg and Alfredsson1998) predict a relaxation length of $0.36a\,Re$. The relaxation process leads to comparatively intense streamwise gradients in the vicinity of the nozzle exit, $x\lesssim 0.01a\,Re$ (Boronin et al. Reference Boronin, Healey and Sazhin2013). Downstream of this region, the streamwise gradient becomes increasingly negligible relative to the radial gradient. The linear analyses in this work thus consider a parallel base flow that is a function of the radial coordinate only, and describe a steady, laminar and axisymmetric liquid jet taken from an analytical solution, previously introduced as profile $U_3$ in Boronin et al. (Reference Boronin, Healey and Sazhin2013). The profile is representative of a nozzle that is sufficiently long to allow the development of the flow into a Hagen–Poiseuille profile.

Throughout this work, $x$, $r$ and $\theta$ denote the axial, radial and azimuthal directions, respectively, and $u$, $v$ and $w$ denote the velocity disturbances in these dimensions. All quantities are normalized by the the radius of the liquid jet, $R$, and its centreline velocity, $U_c$. We further introduce a decomposition of the velocity vector $\boldsymbol {u}=(u,v,w)^\textrm {{T}}$ into a mean, or base flow and a small perturbation component:

(2.1)\begin{equation} \boldsymbol{u}\left(x,r,\theta,t\right)=\bar{\boldsymbol{u}}\left(r\right)+ \boldsymbol{u}'\left(x,r,\theta,t\right). \end{equation}

The analytical expressions for the base flow velocity profiles in the liquid and gas phases, respectively, are

(2.2)\begin{equation} \left.\begin{gathered} \bar{u}_l(r)=\alpha_l \exp \left( \frac{-r^2}{\delta_l^2} \right) + \beta_l, \quad 0 < r < 1,\\ \bar{u}_g(r)=\alpha_g \exp \left\{ \frac{-(r-r_g)^2}{\delta_g^2} \right\}, \quad r > 1. \end{gathered}\right\}\end{equation}

Here and throughout the paper, the indices $g$ and $l$ denote quantities taken from the gas and liquid phases, respectively. Following Boronin et al. (Reference Boronin, Healey and Sazhin2013), we choose $\alpha _l=0.6$, $\beta _l=0.4$, $r_l=0$, $\delta _l=0.7$ and $\delta _g=0.25$. The parameters $\alpha _g$ and $r_g$ are determined so as to satisfy the interfacial continuity of velocity and stress in the streamwise velocity component. We note that the streamwise expansion of the jet balances the radial diffusion. However, the resulting streamwise gradients are small and are thus disregarded. A visualization of the streamwise invariant base flow profile is presented in figure 2.

Figure 2. Base flow of liquid jet in stationary gas. The material interface is located at $r=1$.

The radial extent of the computational domain is three jet radii, $L_r=3$.

We further introduce a set of non-dimensional parameters, including the Reynolds number

(2.3)\begin{equation} Re= \frac{\rho_l U_c R}{\mu_l}, \end{equation}

with $\rho _l$ the fluid density. The Weber number is defined as

(2.4)\begin{equation} We=\frac{\rho_l U_c^2 R}{\sigma}, \end{equation}

with $\sigma$ denoting the interfacial surface tension. The set of non-dimensional parameters describing the two-fluid problem is completed by the ratio of the densities of the two media,

(2.5)\begin{equation} \eta=\frac{\rho_g}{\rho_l}, \end{equation}

and the ratio of their viscosities,

(2.6)\begin{equation} m=\frac{\mu_g}{\mu_l}. \end{equation}

2.2. Linear stability equations

We consider the dynamics of small disturbances to the jet, as described by the linearized Navier–Stokes equations in cylindrical coordinates:

(2.7)\begin{gather} \frac{\partial u'}{\partial t}+\bar{u}\frac{\partial u'}{\partial x} + v' \frac{\partial \bar{u}}{\partial r}= -\frac{\partial p'}{\partial x}+\frac{1}{Re}\nabla^2u' , \end{gather}

(2.8)\begin{gather}\frac{\partial v'}{\partial t}+\bar{u}\frac{\partial v'}{\partial x}= -\frac{\partial p'}{\partial r}+\frac{1}{Re}\left(\nabla^2v' -\frac{v'}{r^2}-\frac{2}{r^2}\frac{\partial w'}{\partial \theta} \right), \end{gather}

(2.9)\begin{gather}\frac{\partial w'}{\partial t}+\bar{u}\frac{\partial w'}{\partial x}= -\frac{1}{r}\frac{\partial p'}{\partial \theta}+\frac{1}{Re}\left( \nabla^2w' + \frac{2}{r^2}\frac{\partial v'}{\partial \theta} - \frac{w'}{r^2} \right) , \end{gather}

(2.10)\begin{gather}\frac{\partial u'}{\partial x} + \frac{1}{r}\frac{\partial}{\partial r}(rv')+\frac{1}{r}\frac{\partial w'}{\partial \theta}=0 , \end{gather}

with the Laplace operator in cylindrical coordinates defined as $\nabla ^2=({\partial ^2}/{\partial x^2})+({1}/{r})(\partial /\partial {r})(r(\partial /\partial {r}))+ ({1}/{r^2})(\partial ^2/\partial \theta ^{2})$. We note that in contrast to formulations in terms of the radial velocity and vorticity, as proposed for instance in Schmid & Henningson (Reference Schmid and Henningson1994), the present formulation, (2.7)–(2.10), is limited to spatial derivatives of second order.

The streamwise-invariant, axisymmetric and steady base flow introduced in § 2.1 allows a normal-mode ansatz for the disturbances, which are assumed to be of the form

(2.11)\begin{equation} \left[u', v', w', p' \right]=\left[ \hat{u}(r), \hat{v}(r), \hat{w}(r), \hat{p}(r) \right] \exp({{\rm i}(\alpha x + n\theta - \omega t)}). \end{equation}

In the same manner, the disturbance of the interface between the two phases $f'$ is expressed as

(2.12)\begin{equation} f'=\hat{f} \exp({{\rm i}(\alpha x + n\theta - \omega t)}). \end{equation}

Here, $\alpha$ and $n$ denote the streamwise and discrete azimuthal wavenumbers, respectively, and $\omega$ represents the frequency of the disturbance. The spatial discretization in the inhomogeneous radial dimension is based on two sets of Chebyshev polynomials that represent the flow in the two phases. The liquid phase uses $N_l=120$ Chebyshev collocation points, and $N_g=70$ collocation points are employed in the gas phase. Substitution of the ansatz (2.11) into the governing equations (2.7)–(2.10) leads to

(2.13)\begin{gather} \left\{ -\frac{1}{Re} \frac{\partial^2}{\partial r^2} -\frac{1}{Re} \frac{1}{r} \frac{\partial }{\partial r} + \left[ -{\rm i}\omega + {\rm i}\alpha \bar{u} +\frac{1}{Re}\left( \frac{n^2}{r^2} + \alpha^2\right) \right] \right\}\hat{u} + {\rm i} \frac{\partial \bar{u}}{\partial r} \hat{v} + {\rm i} \alpha \hat{p}=0, \end{gather}

(2.14)\begin{gather}\left\{ -\frac{1}{Re}\frac{\partial^2 }{\partial r ^2} -\frac{1}{Re}\frac{1}{r}\frac{\partial }{\partial r} + \left[-{\rm i}\omega + {\rm i}\alpha \bar{u} +\frac{1}{Re}\left( \frac{n^2+1}{r^2} + \alpha^2\right) \right] \right\}\hat{v} + \frac{{\rm i}}{Re}\frac{2n}{r^2} \hat{w} + \frac{\partial }{\partial r} \hat{p}=0, \end{gather}

(2.15)\begin{gather}\left\{ -\frac{1}{Re}\frac{\partial^2 }{\partial r ^2} -\frac{1}{Re}\frac{1}{r}\frac{\partial }{\partial r} + \left[ -{\rm i}\omega + {\rm i}\alpha \bar{u} +\frac{1}{Re}\left( \frac{n^2+1}{r^2} + \alpha^2\right) \right] \right\}\hat{w} - \frac{1}{Re}\frac{2n}{r^2} \hat{v} + {\rm i}\frac{n}{r} \hat{p}=0, \end{gather}

(2.16)\begin{gather}{\rm i} \alpha \hat{u} + \frac{1}{r} \frac{\partial }{\partial r} \hat{v} + {\rm i}\frac{n}{r}\hat{w}=0. \end{gather}

The interface disturbance is assumed to be small and is related to the radial disturbance velocity component through the kinematic condition

(2.17)\begin{equation} \frac{{\rm D}f'}{{\rm D}t}=\left( -{\rm i} \omega + {\rm i} \alpha \bar{u}_{l}^I \right) \hat{f}\, \exp({{\rm i}(\alpha x + n\theta - \omega t)}) =\hat{v}_{l}^I \exp({{\rm i}(\alpha x + n\theta - \omega t)}), \end{equation}

where $({{\rm D}}/{{\rm D}t})=(\partial /\partial t)+(\bar {u}(\partial /\partial x))$ denotes the material derivative and the superscript $I$ denotes quantities evaluated at the interface.

2.3. Boundary and interface conditions

The solution of the governing equations (2.13)–(2.16) requires the prescription of boundary and interface conditions. At the outer boundary of the computational domain, in the far field of the gas phase, no-slip and no-penetration boundary conditions are imposed on the eigenfunctions:

(2.18)\begin{equation} \hat{u}_g=0, \quad \hat{v}_g=0,\quad \hat{w}_g=0\quad {\rm at} \ r=L_r. \end{equation}

Special care has to be taken at the centreline because of the singular nature of the cylindrical coordinate system. The assumption of smooth and bounded behaviour of all physical quantities at the centreline, $r \rightarrow 0$, implies that (see e.g. Khorrami, Malik & Ash Reference Khorrami, Malik and Ash1989)

(2.19)\begin{equation} \lim_{r \to 0} \frac{\partial \hat{{\boldsymbol{u}}}}{\partial \theta}=\frac{\partial}{\partial \theta} \left( \hat{u} {\boldsymbol{e}}_z+ \hat{v} {\boldsymbol{e}}_r + \hat{w} {\boldsymbol{e}}_\theta \right)= 0, \end{equation}

and for the pressure

(2.20)\begin{equation} \lim_{r \to 0} \frac{\partial \hat{p}}{\partial \theta}=0. \end{equation}

Substitution of the normal-mode ansatz (2.11) into (2.19)–(2.20) yields

(2.21)\begin{gather} n\hat{v} + \hat{w} = 0, \end{gather}

(2.22)\begin{gather}\hat{v} + n\hat{w} = 0, \end{gather}

(2.23)\begin{gather}n \hat{u} = 0, \end{gather}

(2.24)\begin{gather}n \hat{p}=0. \end{gather}

The resulting centreline conditions depend on the discrete azimuthal number, $n$, as follows:

(2.25)\begin{gather} \frac{\partial\hat{u}}{\partial r} = 0, \quad\hat{v} = 0, \quad \hat{w} = 0, \quad \frac{\partial\hat{p}}{\partial r} = 0 \quad {\rm for} \ n=0, \end{gather}

(2.26)\begin{gather}\hat{u} = 0, \quad -{\rm i} \hat{v} + n \hat{w} = 0, \quad -2 {\rm i} \frac{\partial\hat{v}}{\partial r} + n \frac{\partial\hat{w}}{\partial r} = 0, \quad \hat{p} = 0 \quad {\rm for} \ n=-1 \ {\rm or} \ 1, \end{gather}

(2.27)\begin{gather}\hat{u} = 0, \quad\hat{v} = 0, \quad \hat{w} = 0, \quad \hat{p} = 0 \quad {\rm otherwise}. \end{gather}

We note that because of the linear dependence of conditions (2.21) and (2.22) in the case $|n| = 1$, one of the conditions for that case was derived by applying the continuity equation (2.16).

At the interface ($r = 1$), both the velocity and stress are continuous in the radial dimension. In the parallel velocities and tangential stresses, the differences across the interface are proportional to the surface tension and mean gradient, respectively:

(2.8a–c)\begin{gather} \mathcal{J}(\hat{u})=\mathcal{J}\left( \frac{\partial \bar{u}}{\partial r}\right) \hat{f}, \quad \mathcal{J}(\hat{v})=0, \quad \mathcal{J}(\hat{w})=0, \end{gather}

(2.9a–c)\begin{gather} \mathcal{J}(\widehat\tau_{zr})=\mathcal{J}\left(\mu \frac{\partial^2 \bar{u}}{\partial r^2} \right)\hat{f}, \quad \mathcal{J}(\widehat\tau_{\theta r})=0, \quad \mathcal{J}(\widehat\tau_{rr})=\sigma\nabla^2 \hat{f}.\end{gather}

Here, $\mathcal {J}(\mathcal {X})$ denotes the difference in $\mathcal {X}$ between the two phases at the interface. We note that the relations for $\hat {u}$ and $\hat {\tau }_{zr}$ have been derived by employing first-order Taylor expansions, implying that $\bar {u}_j(r=1+f)=\bar {u}_j^I+({{\rm d}\bar {u}_j^I}/{{\rm d}r}) \hat {f}$. The resulting interface conditions are

(2.30a)\begin{gather} \mathcal{J}\left( {\rm i} \mu \alpha \hat{v} + \mu \frac{\partial \hat{u}}{\partial r} \right)=\mathcal{J}\left( \mu \frac{\partial^2 \bar{u}}{\partial r^2} \right) \hat{f}, \end{gather}

(2.30b)\begin{gather}\mathcal{J}\left( \mu \left(\frac{\partial}{\partial r} - 1 \right)\hat{w} + {\rm i} \mu n\hat{v} \right)=0, \end{gather}

(2.30c)\begin{gather}\mathcal{J}\left( - \frac{\rho}{\rho_1} \hat{p} + \frac{2}{Re}\frac{\mu}{\mu_1} \frac{\partial\hat{v}}{\partial r} \right)=-\frac{1}{We} \left( 1 - n^2 - \alpha^2 \right) \hat{f}. \end{gather}

2.4. Spatial linear analysis

The linearized equations (2.13)–(2.16) permit solutions of both temporally and spatially developing solutions. In this work, we focus exclusively on the spatial case which more faithfully represents the amplification of perturbations in the shear layer of the liquid jet. Spatial stability analysis prescribes a real-valued frequency $\omega$ and azimuthal wavenumber $n$ of the disturbances. As a consequence, only certain complex wavenumbers $\alpha$ satisfy the stability equations (2.13)–(2.16), which gives rise to an eigenvalue problem of the form

(2.31)\begin{equation} {\rm i} \alpha {\boldsymbol{\mathsf{M}}} \hat{\boldsymbol{q}}={\boldsymbol{\mathsf{L}}} \hat{\boldsymbol{q}}. \end{equation}

Here, $\hat {\boldsymbol {q}}$ is the disturbance state vector

(2.32)\begin{equation} \hat{\boldsymbol{q}}=\left[\hat{\phi}_g(r)\enspace \hat{f} \enspace \alpha \hat{f} \enspace \hat{\phi}_l(r)\right]^\textrm{{T}}, \end{equation}

with

(2.33)\begin{equation} \hat{\phi}_i=\left[ \hat{u}_i(r)\enspace \hat{v}_i(r)\enspace \hat{w}_i(r)\enspace \hat{p}_i(r)\enspace \alpha \hat{u}_i(r)\enspace \alpha \hat{v}_i(r)\enspace \alpha \hat{w}_i(r)\right]^\textrm{{T}}. \end{equation}

The operators ${\boldsymbol{\mathsf{M}}}$ and ${\boldsymbol{\mathsf{L}}}$ are given in appendix A.

Since our focus is on the study of disturbances in the absence of exponential instability, we restrict our analyses to exponentially stable settings in which all eigensolutions of the spatial stability problem (2.31) decay in $x$. While the sign of the eigenvalue reveals the stability of the eigensolution in temporal analysis, the nonlinear spatial eigenvalue problem prevents such immediate classification. Existing approaches to the identification of the stability properties of spatial modes involve iterative schemes as applied for instance by Lin & Chen (Reference Lin and Chen1998) and Boronin et al. (Reference Boronin, Healey and Sazhin2013). In the present work, we directly determine the spatial stability of individual eigenfunctions by evaluating their group velocity,

(2.34)\begin{equation} c_g\equiv \frac{\partial \omega}{\partial \alpha}, \end{equation}

using information provided by the adjoint stability equations. More specifically, a positive spatial growth rate, i.e. a negative imaginary part of $\alpha$, implies spatial growth if the mode travels downstream, indicated by a positive group velocity. If the group velocity of the eigenfunction is negative, it travels upstream and thus undergoes decay. Analogously, a mode with a negative spatial growth rate, i.e. a positive imaginary part of $\alpha$, will decay if its group velocity is positive, and amplify if its group velocity is negative.

In this context, we define the inner product as the volume integrated dot product within the computational domain:

(2.35)\begin{equation} \langle \boldsymbol{a}, \boldsymbol{b}\rangle=\int_\varOmega \boldsymbol{a}^{\mathsf{H}} \boldsymbol{b} \, {\rm d} V, \end{equation}

with superscript ${\mathsf{H}}$ denoting the complex conjugate transpose. The adjoint eigenfunctions, $\hat {\boldsymbol {q}}^{{\dagger}ger }$, are the solution of the adjoint eigenvalue problem

(2.36)\begin{equation} \alpha^*{\boldsymbol{\mathsf{M}}}^{\mathsf{H}}\hat{\boldsymbol{q}}^{{\dagger}ger}={\boldsymbol{\mathsf{L}}}^{\mathsf{H}} \hat{\boldsymbol{q}}^{{\dagger}ger}, \end{equation}

where superscript $*$ denotes complex conjugate. Solutions to the forward and adjoint eigenvalue problem satisfy a bi-orthogonality condition, so that upon proper normalization

(2.37)\begin{equation} \langle \hat{\boldsymbol{q}}^{{\dagger}ger}_i, {\boldsymbol{\mathsf{M}}} \hat{\boldsymbol{q}}_j \rangle=\delta_{ij}. \end{equation}

We derive a relation for the group velocity by applying the partial derivative with respect to $\omega$ to (2.31), and integrating over the computational domain:

(2.38)\begin{equation} \frac{\partial \alpha}{\partial \omega} \langle \hat{\boldsymbol{q}}^{{\dagger}ger}, {\boldsymbol{\mathsf{M}}}\hat{\boldsymbol{q}} \rangle=\left\langle \hat{\boldsymbol{q}}^{{\dagger}ger}, \frac{\partial {\boldsymbol{\mathsf{L}}}}{\partial \omega}\hat{\boldsymbol{q}} \right\rangle. \end{equation}

Application of the bi-orthogonality condition (2.37) to (2.38) leads to

(2.39)\begin{equation} \frac{\partial \alpha}{\partial \omega}=\left\langle \hat{\boldsymbol{q}}^{{\dagger}ger}, \frac{\partial {\boldsymbol{\mathsf{L}}}}{\partial \omega}\hat{\boldsymbol{q}} \right\rangle. \end{equation}

The group velocity of an eigenfunction $\hat {\boldsymbol {q}}$, and by extension its exponential stability within the spatial stability framework, may thus be evaluated by the inner product with its corresponding adjoint eigenvector $\hat {\boldsymbol {q}}^{{\dagger}ger }$:

(2.40)\begin{equation} \frac{\partial \omega}{\partial \alpha} =\left\langle \hat{\boldsymbol{q}}^{{\dagger}ger}, \frac{\partial {\boldsymbol{\mathsf{L}}}}{\partial \omega}\hat{\boldsymbol{q}}\right\rangle^{-1}. \end{equation}

We note that while iterative approaches require multiple solutions of the stability problem for each individual mode, the gradient information provided by the adjoint problem identifies the group velocity of all eigensolutions after solving the adjoint eigenvalue problem (2.36) once.

2.5. Optimization problem

Even when all the eigensolutions of the spatial stability problem (2.31) are stable, the non-normal nature of the linearized Navier–Stokes equations provides a pathway for disturbance growth. The study of this non-modal, or transient, amplification of disturbances can be facilitated by considering an optimization problem in a suitably chosen norm. Earlier analyses of non-modal growth in the two-fluid setting followed studies in single-fluid shear flows in effectively optimizing the kinetic energy of disturbances (e.g. Renardy Reference Renardy1987; Yecko & Zaleski Reference Yecko and Zaleski2005; Boronin et al. Reference Boronin, Healey and Sazhin2013). In the wider context of flow atomization, the amplification of distortions specifically of the material interface is nonetheless of greater relevance.

We investigate the specific conditions that amplify distortions of the interface by introducing an interface energy norm of the form

(2.41)\begin{equation} ||\boldsymbol{q}||_{IE}=\int_{0}^{{2{\rm \pi}}/{n}} \int_{0}^{{2{\rm \pi}}/{\omega}} \frac{1}{We} (-1+n^2+\alpha^2) \hat{f}^* \hat{f} \, {\rm d}t \, {\rm d}\theta. \end{equation}

The norm represents a measure for the surface tension energy stored in the material interface owing to the perturbation displacement $\hat {f}$. The derivation is presented in appendix B. Together with the norm (2.41), we introduce the suitably chosen weight matrix ${\boldsymbol{\mathsf{F}}}$ which allows the transformation of the interface energy norm into a 2-norm so that $||\boldsymbol {q}||_{IE}\equiv ||{\boldsymbol{\mathsf{F}}}\boldsymbol {q}||_2$. We note that the interface energy norm (2.41) partitions the state vector and as such describes a non-trivial semi-norm. As a consequence, ${\boldsymbol{\mathsf{F}}}$ is rank-deficient. Specifically, the velocity and pressure components of both the gas and liquid phases are within the kernel of ${\boldsymbol{\mathsf{F}}}$. The formulation of a well-posed optimization problem thus calls for the introduction of a full norm (see e.g. Foures, Caulfield & Schmid Reference Foures, Caulfield and Schmid2012), which we choose as the total energy norm:

(2.42)\begin{align} ||\boldsymbol{q}||_E &= \frac{1}{2} \eta \int_{1}^{L_r} \int_{0}^{{2{\rm \pi}}/{n}} \int_{0}^{{2{\rm \pi}}/{\omega}} \left(\hat{u}_g^* \hat{u}_g + \hat{v}_g^* \hat{v}_g+ \hat{w}_g^* \hat{w}_g\right)r\,{\rm d}t\,{\rm d}\theta\,{\rm d}r\nonumber\\ &\quad + \frac{1}{2} \int_{0}^{1} \int_{0}^{{2{\rm \pi}}/{n}} \int_{0}^{{2{\rm \pi}}/{\omega}} \left( \hat{u}_l^* \hat{u}_l + \hat{v}_l^* \hat{v}_l + \hat{w}_l^* \hat{w}_l \right)r\,{\rm d}t\,{\rm d}\theta\,{\rm d}r \nonumber\\ &\quad + \int_{0}^{{2{\rm \pi}}/{n}} \int_{0}^{{2{\rm \pi}}/{\omega}} \frac{1}{We} (-1 + n^2+\alpha^2) \hat{f}^* \hat{f}\, {\rm d}t \, {\rm d}\theta. \end{align}

Analogous to the interface energy norm, the energy norm can be converted to a 2-norm by means of the matrix ${\boldsymbol{\mathsf{F}}}_N$, with $||\boldsymbol {q}||_E\equiv ||{\boldsymbol{\mathsf{F}}}_N\boldsymbol {q}||_2$. Since (2.42) covers all components of the state vector, the null space of ${\boldsymbol{\mathsf{F}}}_N$ is empty.

Our optimization problem seeks to find the specific, optimal, initial condition $\boldsymbol {q}_0$, which over a prescribed spatial distance evolves into the particular final condition $\boldsymbol {q}_1$ that induces the most energetic distortions of the jet interface, as measured by the norm $\|\boldsymbol {q}\|_{IE}$. In our solution of the optimization problem, we follow Hack & Moin (Reference Hack and Moin2017) in expressing the initial condition as a weighted sum of eigenfunctions of the two-fluid stability problem (2.31):

(2.43)\begin{equation} \boldsymbol{q}_0={\boldsymbol{\mathsf{Q}}}_0 \boldsymbol{\kappa}. \end{equation}

Here, the columns of ${\boldsymbol{\mathsf{Q}}}_0 \in \mathbb {C}^{N_Q\times N}$, where $N_Q=7N_{l}+7N_{g}+2$, contain the $N$ least stable eigenfunctions $\hat {\boldsymbol {q}}_{j,0}$ and $\boldsymbol {\kappa }\in \mathbb {C}^{N\times 1}$ assigns each of the modes a specific weight. Consistently, the solution at the final position $x_1$ may be expressed as

(2.44)\begin{equation} \boldsymbol{q}_1={\boldsymbol{\mathsf{Q}}}_1\boldsymbol{\kappa}, \end{equation}

with the columns of ${\boldsymbol{\mathsf{Q}}}_1$ representing the evolution of the $n$ eigenfunctions:

(2.45)\begin{equation} \hat{\boldsymbol{q}}_{j,1}=\hat{\boldsymbol{q}}_{j,0}\exp\left({\rm i}\alpha_j \left(x_1-x_0\right)\right). \end{equation}

The linear nature of the problem implies that $\boldsymbol {\kappa }$ is identical in both expressions. Since ${\boldsymbol{\mathsf{Q}}}_0$ has full column rank, we can explicitly determine its Moore–Penrose pseudo-inverse ${\boldsymbol{\mathsf{Q}}}_0^+=({\boldsymbol{\mathsf{Q}}}_0^{{\mathsf{H}}}{\boldsymbol{\mathsf{Q}}}_0)^{-1}{\boldsymbol{\mathsf{Q}}}_0^{{\mathsf{H}}}$, leading to

(2.46)\begin{equation} {\boldsymbol{\mathsf{Q}}}_0^+ \boldsymbol{q}_0=\boldsymbol{\kappa}. \end{equation}

Substitution into (2.44) yields

(2.47)\begin{equation} \boldsymbol{q}_1= {\boldsymbol{\mathsf{L}}}_1 \boldsymbol{q}_0, \end{equation}

where ${\boldsymbol{\mathsf{L}}}_1={\boldsymbol{\mathsf{Q}}}_1 {\boldsymbol{\mathsf{Q}}}_0^+$ is a pseudo-propagator that advances arbitrary initial disturbances $\boldsymbol {q}_0$ within the column space of ${\boldsymbol{\mathsf{Q}}}_0$ from $x_0$ to $x_1$. Without loss of generality, we choose $x_0=0$ in this work.

The objective of finding the specific initial condition $\boldsymbol {q}_0$ that maximizes the amplification of interface distortions with frequency $\omega$ at a given target location $x_1$ may thus be expressed as the functional

(2.48)\begin{align} G(\omega,x_1)&= \frac{1}{2} \max_{\boldsymbol{q}_0} \frac{||\boldsymbol{q}_1||_{IE}^2}{||\boldsymbol{q}_0||_E^2} \end{align}

(2.49)\begin{align} &= \frac{1}{2} \max_{\boldsymbol{q}_0} \frac{||{\boldsymbol{\mathsf{L}}}_1 \boldsymbol{q}_0||_{IE}^2}{||{\boldsymbol{\mathsf{L}}}_0 \boldsymbol{q}_0||_E^2} \end{align}

(2.50)\begin{align} &= \frac{1}{2} ||{\boldsymbol{\mathsf{F}}}{\boldsymbol{\mathsf{L}}}_1 ({\boldsymbol{\mathsf{F}}}_N{\boldsymbol{\mathsf{L}}}_0)^+||_2^2. \end{align}

Here, $\boldsymbol {q}_1=\boldsymbol {q}(x_1)$ and the introduction of ${\boldsymbol{\mathsf{L}}}_0={\boldsymbol{\mathsf{Q}}}_0{\boldsymbol{\mathsf{Q}}}_0^+$ into the denominator restricts $\boldsymbol {q}_0$ to the column space of ${\boldsymbol{\mathsf{Q}}}_0$. We further note that the majority of the results presented in the following focus on the maximum gain, which represents the maximum $G$ for all possible perturbation frequencies and target locations, $G_{max}=\max _{\omega ,x_1}G(\omega ,x_1)$.

The objective functional $G$ relates the interface energy at the target position $x_1$ to the sum of the kinetic and interface perturbation energies at the initial position $x_0$. In other words, $G$ quantifies the interface energy at the target location normalized by the total initial perturbation energy. In doing so, it prescribes a fair and conservative measure for the amplification of interface distortions as $G$ attains values larger than unity only if the perturbation interface energy at the considered target location exceeds the total energy at the initial location $x_0$, i.e. if energy has been extracted from the mean shear and redistributed into the potential energy of the interface. A trivial redistribution of initial perturbation kinetic energy into perturbation interface energy at the target location therefore does not add to the gain. We also note that the pre-factor of $\frac {1}{2}$ has been added for convenience only. In its absence, the critical value separating genuine amplification from decay or trivial redistribution would be $2$. In what follows, we apply the conceptual framework presented above in the demonstration of the linear amplification of distortions to the interface of a liquid jet by a non-exponential mechanism.

3. Mechanism for the amplification of interface distortions

In the following, we demonstrate a mechanism that allows the amplification of distortions to the material interface of a liquid jet in exponentially stable settings. In a first step, we analyse the linear stability of a jet described by profile (2.2) by solving the spatial stability problem (2.31) for the considered parameters $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $n=1$ and $\omega =0.7$. We note that the considered viscosity and density ratios are comparable to those of the combination of water and air. Figure 3(a) presents the real and imaginary parts of the computed eigenvalues in terms of the complex streamwise wavenumber, $\alpha$, of the disturbances. The polynomial nature of the spatial stability problem impedes a classification of the exponential stability of eigensolutions based on the sign of the eigenvalue alone, as discussed in § 2.4. The determination of the stability of a given eigensolution instead has to consider both the eigenvalue and the associated group velocity. For the present setting, figure 3(b) presents the group velocities of the individual eigenvalues along with the imaginary parts of the wavenumber. The group velocity is positive for all considered eigenvalues, and the two eigensolutions with negative $\alpha _i$ are thus exponentially unstable while all other modes with $\alpha _i>0$ are stable.

Figure 3. (a) Spectrum of the complex eigenvalue $\alpha$ for $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $n=1$ and $\omega =0.7$. (b) Group velocity $c_g$ and growth rate $\alpha _i$ of the eigenfunctions. Blue and red circles mark the unstable modes.

The imaginary parts of the eigenvalues, $\alpha _i$, of the respective least stable modes at eight different azimuthal numbers, $n$, are presented in figure 4 as a function of the disturbance frequency $\omega$. The parameters are $Re=7500$, $We=3000$, $m=0.01$ and $\eta =0.001$, and the considered azimuthal wavenumbers range from $n=0$ to $n=7$. In all cases, the corresponding group velocities are positive, and the sign of the eigenvalues thus determines the stability of the solutions. The modes computed for $n=0$ to $n=4$ are exponentially unstable for certain perturbation frequencies. The highest amplification occurs for $n=1$, which is referred to as the sinuous mode of the K–H instability (Söderberg Reference Söderberg2003). For $n > 4$, the most unstable modes decay for all considered frequencies. Our analyses of the transient amplification of interface perturbations thus focus on this wavenumber range. Analogous to the optimal modes with non-zero transverse wavenumber considered by Yecko & Zaleski (Reference Yecko and Zaleski2005), these non-axisymmetric modes directly provide a pathway for the generation of transverse perturbations.

Figure 4. Spatial growth rate $\alpha _i$ as a function of the disturbance frequency $\omega$ for the least stable modes at azimuthal wavenumbers $n=0$ (blue), $n=1$ (red), $n=2$ (yellow), $n=3$ (purple), $n=4$ (green), $n=5$ (black solid), $n=6$ (black dashed) and $n=7$ (black dotted). The parameters are $Re=7500$, $We=3000$, $m=0.01$ and $\eta =0.001$.

The solution of the optimization problem in the interface energy norm $\|\boldsymbol {q}\|_{IE}$ establishes that a significant amplification of interfacial distortions can occur in the absence of exponential instability. The 2-norm of the operator ${\boldsymbol{\mathsf{F}}}{\boldsymbol{\mathsf{L}}}_1({\boldsymbol{\mathsf{F}}}_N{\boldsymbol{\mathsf{L}}}_0)^+$ is determined by the singular value decomposition:

(3.1)\begin{equation} {\boldsymbol{\mathsf{U}}}{\boldsymbol{\varSigma}}{\boldsymbol{\mathsf{V}}}^{\mathsf{H}} = {\boldsymbol{\mathsf{F}}}{\boldsymbol{\mathsf{L}}}_1({\boldsymbol{\mathsf{F}}}_N{\boldsymbol{\mathsf{L}}}_0)^+. \end{equation}

Here, ${\boldsymbol{\mathsf{U}}}$ is an orthogonal matrix whose columns are the left singular vectors, and the columns of ${\boldsymbol{\mathsf{V}}}$ contain the right singular values. The entries of the real, positive definite diagonal matrix ${\boldsymbol {\varSigma }}$ describe the singular values $\sigma _j$. Specifically, the largest singular value describes the gain, $G\equiv \sigma _{max}$, and the associated right singular vector gives the corresponding initial state. The parameters are $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $n=5$ and $\omega =2.2$ and the target location is $x_1=11.8$. The amplification of the distortion to the material interface is quantified by evaluating the magnitude of the interface displacement, $f$. The absolute value of the interface distortion, normalized by its initial value, is presented in figure 5(a) as a function of the streamwise coordinate. The maximum amplification ratio of the interface distortion at $x=11.8$ is $|\,f(x)/\,f_0|\simeq 1042$. The growth of the interface distortion towards the target location, $x_1$, shows an oscillatory, non-monotonic behaviour. The streamwise evolution of the integrated kinetic energy in the velocity perturbations and surface potential energy is presented in figure 5(b,c). The kinetic energy in all components generally increases from the initial location $x_0$ to the target location $x_1=11.8$. While the energy in the streamwise components gains most, we note that the azimuthal and critically the radial components also amplify. In contrast to the streamwise and azimuthal components, the gain in the radial perturbations is non-monotonic, and the energy in that component is below its peak at the target location. Analogously, the trend in the interface potential energy is non-monotonic, although the peak of this quantity coincides with the target location. As we demonstrate in § 4, the non-monotonic increase in both the interface displacement and the kinetic energy in the radial velocity perturbations is related to a temporary redistribution between the potential energy in the interface due to surface tension and the kinetic energy of the velocity perturbations.

Figure 5. (a) Spatial evolution of the absolute value of the ratio of interface disturbance amplitude. (b) Spatial evolution of the integrated kinetic energy in the streamwise velocity component (blue), radial velocity component (red) and azimuthal component (yellow), their sum (black) and surface tension energy of the interface perturbation (green). (c) Detail showing the spatial evolution of the integrated kinetic energy in the radial velocity component (red) and surface tension energy of the interface perturbation (green). The dashed lines mark the target location, $x_1=11.8$. The parameters are $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $n=5$ and $\omega =2.2$.

The significance of the non-exponential amplification of the interface distortion is established by comparing with the fastest-growing exponential instability, which occurs for an azimuthal wavenumber $n=1$ (see figure 4). The growth rate of the mode is $\alpha _i=0.1028$, leading to an expected growth of the most unstable mode during the considered interval of $\exp (\alpha _i x_1)=\exp (0.1028 \times 11.8)=3.364$, which is significantly smaller than the ratio of $|\,f(x)/\,f_0|\simeq 1042$ reported for the optimal interface distortion. This result demonstrates that over finite streamwise intervals, the algebraic amplification of interface distortions by means of the presented mechanism may considerably exceed that by the exponential K–H instability.

A visualization of the spatial evolution of the disturbance field is presented in figure 6. Vectors mark the streamwise and radial perturbation velocities, $u'$ and $v'$, in the $x$–$r$ plane and colour contours show the azimuthal component, $w'$, of the perturbation velocity. The growth in the intensity of the perturbation velocities coincides with a reorientation of the initially tilted fluid structure by the mean shear. As a consequence, the distortion of the material interface is amplified and peaks at $x=11.8$, marked by the vertical dashed line. Both the tilting of the perturbation structures and the amplification of velocity perturbations normal to the direction of the mean shear are characteristic of the Orr mechanism, as sketched in figure 1. In the two-fluid setting, these perturbations are also normal to the the material interface, and thus induce an amplification of interface distortions.

Figure 6. Side view of the streamwise evolution of the optimal initial condition. Vectors of the streamwise and radial disturbance fields, $u'$ and $v'$, with colour contours of the azimuthal disturbance field, $w'$. The red line marks the position of material interface. The white dashed line marks the target location, $x_1$. The parameters are $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $n=5$ and $\omega =2.2$.

We note that the presented outcome qualitatively differs from that of the lift-up mechanism which generates low-frequency perturbations in the streamwise perturbation velocity component (e.g. Gustavsson Reference Gustavsson1991; Hack & Zaki Reference Hack and Zaki2015) and which was identified in earlier analyses of non-exponentially amplifying disturbances in two-fluid settings (Yecko & Zaleski Reference Yecko and Zaleski2005; Boronin et al. Reference Boronin, Healey and Sazhin2013). Owing to the kinematic coupling of the interface distortion to the velocity component normal to the interface, these perturbations are, however, unlikely to meaningfully advance the atomization process. The present analysis of perturbations which induce a significant amplification of the distortion of the material interface by means of the Orr mechanism is facilitated by the choice of the objective functional in terms of the perturbation interface energy (2.41). This approach differs from the optimization of the total energy of the perturbations, as considered in previous studies (e.g. South & Hooper Reference South and Hooper1999; Yecko & Zaleski Reference Yecko and Zaleski2005).

4. Energy budget

In the following, we seek to provide a more comprehensive documentation of the presented mechanism for amplifying interface distortions by examining the budget of the perturbation kinetic energy. The energy budget is derived by taking the dot product of the linearized governing equations, (2.7)–(2.9), with the velocity perturbation vector for each phase and substituting the continuity equation, (2.10), as well as the interface coupling condition (2.29). Our analyses of the energy budget will focus on integral quantities, defined in terms of the spatio-temporal integral

(4.1)\begin{equation} \int_{\varOmega}\chi\,\textrm{d}\varOmega= \int_{0}^{x}\int_{0}^{{2{\rm \pi}}/{\omega}}\int_{0}^{{2{\rm \pi}}/{n}}\int_{0}^{1} \chi r\,{\rm d}r \, {\rm d}\theta \, {\rm d}t \, {\rm d}x. \end{equation}

Similarly, spatio-temporal integration on the jet interface for a quantity of interest, $\chi$, is denoted as

(4.2)\begin{equation} \int_{\mathcal{S}}\chi\,\textrm{d}{\mathcal{S}}=\int_{0}^{x}\int_{0}^{{2{\rm \pi}}/{\omega}}\int_{0}^{{2{\rm \pi}}/{n}} \left[ \chi \right]_{r=1} \, \textrm{d}\theta \, \textrm{d}t \, \textrm{d}x. \end{equation}

The linearized energy budgets for the three perturbation velocity components thus are

(4.3)\begin{align} \int_{\varOmega}\frac{\textrm{D}}{\textrm{D}t}\left( \frac{1}{2}u'^2 \right)\textrm{d}\varOmega &= \underbrace{\int_{\varOmega}p' \frac{\partial u'}{\partial x}\,\textrm{d}\varOmega}_{\substack{\text{Pressure-velocity} \\ \text{correlation}}} + \underbrace{\int_{\varOmega}\frac{\partial}{\partial x}\left\{ -p'u'+\frac{1}{Re}\frac{\partial }{\partial x} \left(\frac{1}{2}u'^2 \right) \right\}\textrm{d}\varOmega}_\text{Convection} \nonumber\\ &\quad + \underbrace{\int_{\mathcal{S}}\frac{1}{Re} \left( -u' \frac{\partial v'}{\partial x} -u'f' \frac{\partial^2 \bar{u}_l}{\partial r^2} \right)\textrm{d}\mathcal{S}}_\text{Liquid viscous stress} - \underbrace{\int_{\varOmega}u'v'\frac{\partial \bar{u}}{\partial r}\,\textrm{d}\varOmega}_\text{Production} \nonumber\\ &\quad - \underbrace{\int_{\varOmega}\frac{1}{Re} \left( \frac{\partial u'}{\partial x}\frac{\partial u'}{\partial x} + \frac{1}{r}\frac{\partial u'}{\partial \theta} \frac{1}{r}\frac{\partial u'}{\partial \theta} + \frac{\partial u'}{\partial r}\frac{\partial u'}{\partial r} \right)\textrm{d}\varOmega}_\text{Dissipation} \nonumber\\ &\quad + \underbrace{\int_{\mathcal{S}}\frac{1}{Re} m\left( u' \frac{\partial v'_g}{\partial x} + u' \frac{\partial u'_g}{\partial r} + u'f' \frac{\partial^2 \bar{u}_g}{\partial r^2} \right)\textrm{d}{\mathcal{S}}}_\text{Gas viscous stress}, \\[-2.3pc]\nonumber\end{align}

(4.4)\begin{align} &\int_{\varOmega}\frac{\textrm{D}}{\textrm{D}t}\left( \frac{1}{2}v'^2 \right)\textrm{d}\varOmega \nonumber\\ &\quad= \underbrace{\int_{\varOmega}\frac{p'}{r}\frac{\partial }{\partial r}(rv')\,\textrm{d}\varOmega}_{\substack{\text{Pressure-velocity} \\ \text{correlation}}} + \underbrace{\int_{\mathcal{S}} \left( -\eta p'_g v'_g + \frac{2m}{Re}v\frac{\partial v'_g}{\partial r} \right)\textrm{d}{\mathcal{S}}}_\text{Gas viscous stress} - \underbrace{\int_{\mathcal{S}}\frac{1}{Re}v'\frac{\partial v'}{\partial r}\,\textrm{d}\mathcal{S}}_\text{Liquid viscous stress} \nonumber \nonumber\\ &\qquad - \underbrace{\int_{\varOmega}\frac{1}{Re} \left( \frac{\partial v'}{\partial r}\frac{\partial v'}{\partial r} + \frac{1}{r}\frac{\partial v'}{\partial \theta}\frac{1}{r}\frac{\partial v'}{\partial \theta} + \frac{\partial v'}{\partial x}\frac{\partial v'}{\partial x} + \frac{v'^2}{r^2} + \frac{2}{r^2}v'\frac{\partial w'}{\partial \theta} \right)\textrm{d}\varOmega}_\text{Dissipation} \nonumber\\ &\qquad + \underbrace{\int_{\varOmega}\frac{1}{Re} \frac{\partial^2}{\partial x \partial x} \left( \frac{1}{2}v'^2 \right)\textrm{d}{\varOmega{}}}_\text{Convection} - \underbrace{\int_{\mathcal{S}}\frac{1}{We} \left( -1 + n^2 + \alpha^2 \right)f'v'\,\textrm{d}\mathcal{S}}_\text{Surface tension energy}, \end{align}

(4.5)\begin{align} & \int_{\varOmega}\frac{\textrm{D}}{\textrm{D}t}\left( \frac{1}{2}w'^2 \right)\textrm{d}\varOmega \nonumber\\ &\quad = \underbrace{\int_{\varOmega}\frac{p'}{r}\frac{\partial w'}{\partial \theta}\,\textrm{d}\varOmega}_{\substack{\text{Pressure-velocity} \\ \text{correlation}}} + \underbrace{\int_{\varOmega}\frac{1}{Re} \frac{\partial^2}{\partial x \partial x} \left( \frac{1}{2}w'^2 \right) \textrm{d}\varOmega}_\text{Convection} + \underbrace{\int_{\mathcal{S}}\frac{1}{Re} \left( -w'\frac{\partial v'}{\partial \theta} -w'^2 \right)\textrm{d}\mathcal{S}}_\text{Liquid viscous stress} \nonumber \nonumber\\ &\qquad - \underbrace{\int_{\varOmega}\frac{1}{Re} \left( \frac{\partial w'}{\partial r}\frac{\partial w'}{\partial r} + \frac{1}{r}\frac{\partial w'}{\partial \theta}\frac{1}{r}\frac{\partial w'}{\partial \theta} + \frac{\partial w'}{\partial x}\frac{\partial w'}{\partial x} + \frac{w'^2}{r^2} - \frac{2}{r^2}w'\frac{\partial v'}{\partial \theta} \right)\textrm{d}\varOmega}_\text{Dissipation} \nonumber\\ &\qquad + \underbrace{\int_{\mathcal{S}}\frac{1}{Re} m \left( w'\frac{\partial w'_g}{\partial r} + w'\frac{\partial v'_g}{\partial \theta} - w'w'_g \right)\textrm{d}\mathcal{S}}_\text{Gas viscous stress}. \end{align}

Unless marked by a subscript $g$, all quantities are evaluated in the liquid phase. In this context, we also note that the evolution of the surface tension energy is related to the surface tension term in the energy budget (4.6) as follows:

(4.6)\begin{equation} \int_{\mathcal{S}}\frac{\textrm{D}}{\textrm{D}t}\left( \frac{1}{We}(-1+n^2+\alpha^2)f'^2 \right)\textrm{d}\mathcal{S} = 2\int_{\mathcal{S}}\frac{1}{We} \left(-1+n^2+\alpha^2 \right)f'v'\,\textrm{d}\mathcal{S}. \end{equation}

The evolution of the budget terms for the three velocity components from the initial location $x=0$ to the target location $x_1=11.8$ is presented in figure 7. Consistent with the case considered in the previous section, the parameters are $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $n=5$ and $\omega =2.2$. Positive values indicate that the respective term introduces perturbation kinetic energy in the considered component while negative-valued terms extract perturbation kinetic energy from the component. The sum of all quantities at a given downstream location denotes the variation of perturbation kinetic energy in the corresponding component. Consistent with the considered optimization problem, the energy budgets are normalized by the total perturbation energy at the initial location, $x=0$.

Figure 7. Evolution of the terms of the component-wise energy budgets. (a) Streamwise component. (b) Radial component. (c) Azimuthal component. Production (magenta), gas viscous stress (purple), liquid viscous stress (green), PVC (blue), surface tension energy (red), dissipation (orange) and convective (yellow) terms, and their sum (back dotted). The vertical dashed line marks the target location, $x_1=11.8$. The parameters are $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $n=5$ and $\omega =2.2$.

When considering the budget of the streamwise velocity component in figure 7(a), the production term is positive throughout the evolution of perturbation. This term acts as the main energy source of the perturbations by extracting energy from the mean shear. The gas viscous stress terms are of positive sign in both the streamwise and radial dimensions, although their magnitudes are small compared with the production term. In contrast, the dissipation terms are negative for all three components, meaning that energy is drained by viscous dissipation.

Of particular interest is the pressure–velocity correlation (PVC) term. The sum of the PVC terms of the individual velocity components is zero, and the PVC terms thus introduce an inter-component redistribution of perturbation energy. For the present case, the PVC term in the budget of the streamwise velocity component remains negative for the entire process, while positive values are observed for the radial component. In the azimuthal component, the PVC term is initially marginally negative, and eventually attains positive values for $x\gtrsim 4$. The PVC term thus drains energy from the streamwise component and reintroduces it into the radial and azimuthal components. In this manner, perturbation kinetic energy is introduced into the streamwise component by the mean shear, and transferred to the azimuthal and radial components via the PVC terms.

The surface tension energy term only appears in the budget of the radial component, (4.4), and the growth of the energy stored in the interface energy thus occurs independent of the budgets representing the streamwise and azimuthal components. The perturbation energy conveyed to the interface can therefore be studied by considering the interfacial stress condition, (2.30c). Multiplication by $v'$ and integration along the liquid surface yields

(4.7)\begin{align} &\underbrace{\int_{\mathcal{S}}-p'_l v'\,\textrm{d}\mathcal{S}}_{\substack{\text{Liquid} \\ \text{pressure-velocity} \\ \text{correlation}}} + \underbrace{\int_{\mathcal{S}}\frac{2}{Re} \left( v' \frac{\partial v'_l}{\partial r} \right) \textrm{d}\mathcal{S}}_{\substack{\text{Liquid} \\ \text{viscous stress}}} + \underbrace{\int_{\mathcal{S}}\eta p'_g v'\,\textrm{d}\mathcal{S}}_{\substack{\text{Gas} \\ \text{pressure-velocity} \\ \text{correlation}}} - \underbrace{\int_{\mathcal{S}}\frac{2m}{Re} \left( v' \frac{\partial v'_g}{\partial r} \right) \textrm{d}\mathcal{S}}_{\substack{\text{Gas} \\ \text{viscous stress}}} \nonumber\\ &\quad + \underbrace{\int_{\mathcal{S}}\frac{1}{We} \left( -1 + n^2 + \alpha^2 \right)f'v'\,\textrm{d}\mathcal{S}}_{\substack{\text{Surface tension} \\ \text{energy}}} =0. \end{align}

The terms in (4.7) are presented in figure 8. The sum of all terms is zero at all downstream locations, and the gas PVC term is of negligible amplitude. Any gain in the surface perturbation energy can thus be attributed to the balance of the remaining three terms. The results show that the main source of energy is the liquid PVC term. During the early stages ($x\leqslant 7$), the magnitude of the liquid PVC term effectively matches that of the energy absorbed by the surface tension energy while viscous stress terms are nearly zero. Past this point, the liquid PVC and the surface tension terms grow rapidly with matching oscillation patterns. In contrast, the viscous terms grow smoothly, again with matching trends in the two liquids. The positive contribution by the gas viscous stress is moderately exceeded by the negative contribution of the liquid viscous stress, with the difference being compensated by the liquid PVC term. The outcome establishes that the gain in surface tension energy is indeed driven by the liquid PVC term.

Figure 8. Evolution of the terms of the energy budget at the liquid jet interface. Liquid PVC (blue), gas viscous stress (yellow), liquid viscous stress (green), gas PVC (purple) and surface tension energy (red) terms, and their sum (black dotted line). The black dashed line marks the target location, $x_1=11.8$. The parameters are $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $n=5$ and $\omega =2.2$.

In summary, the discussion of the energy budget establishes the mechanics of the algebraic amplification of interface distortions. Perturbation kinetic energy is generated by the mean shear in the streamwise perturbation velocity component and transferred to the radial component through the PVC term. Since only the radial component is kinematically coupled to the interface displacement, this redistribution mechanism is critical in enabling the amplification of interface distortions. At the interface, the energy introduced by the PVC term acts as a source for increasing the surface tension energy.

After the qualitative demonstration of the algebraic amplification of interface distortions, and the analysis of the underlying dynamics, we seek in the following to quantify the scaling of the mechanism with several key parameters.

5. Parameter influence and scaling

The considered problem is governed by five parameters: the Reynolds number $Re$, the Weber number $We$, the ratio of the densities of the two fluids $\eta$, the ratio of their viscosities $m$ and the azimuthal wavenumber of the perturbations $n$. The frequency of the perturbations $\omega$ and the target location $x_1$ of the analysis of optimal disturbances describe two additional parameters in the determination of the gain $G$. As noted in § 2.5, these parameters are nonetheless eliminated in the consideration of the maximum gain $G_{max}$ which is defined as the maximum of all $G(\omega ,x_1)$. Following its definition in (2.50), the objective functional $G$ relates the interface energy at the target location $x_1$ to the sum of the kinetic and interface perturbation energies at the initial position $x=0$. In other words, $G$ quantifies the interface energy at $x_1$ normalized by the total initial perturbation energy. In doing so, it defines a conservative measure for the gain of the interface energy that discounts the trivial redistribution of initial perturbation kinetic energy into surface tension energy at the target location. We note that a side effect of this approach is that the initial value of $G$ at $x=0$ is in general not unity.

The influence of the target location $x_1$ on the objective functional $G$ is presented in figure 9. We note that the presented gain in surface tension energy at each target location is the result of an independent evaluation of the optimization problem. For example, the maximum $G$ for the parameters $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $\omega = 2.2$ and $n = 5$ is $G=2.638$ at $x_1=11.8$ in figure 9.

Figure 9. Energy gain $G$ as a function of the target location $x_1$. The parameters are $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $\omega = 2.2$ and $n = 5$.

In the remainder, we will focus on the maximum gain $G_{max}=\max _{\omega ,x_1}G(\omega , x_1)$ which is obtained by maximizing $G$ for all possible combinations of $\omega$ and $x_1$. The azimuthal wavenumber is chosen from the range of exponentially stable settings, $n=\{5,6,7\}$. We maintain a constant base flow and investigate the effect of the Reynolds and Weber numbers as well as the viscosity and density ratios of the two fluids on the amplification of perturbations.

An overview of the effect of the Reynolds and Weber numbers on $G_{max}$ is presented in figure 10. The remaining parameters are kept constant at $m=0.01, \eta =0.001$ and $n=5$. The effect of the Reynolds number is monotonic so that higher $Re$ also leads to higher gain of the objective functional. On the other hand, an increase of the Weber number leads to a local maximum in $G_{max}$, with further increases of $We$ reducing the gain. More generally speaking, the presented result establishes that both viscosity and surface tension have an appreciable effect on the amplification of interface distortions.

Figure 10. Contours of the maximum energy gain $G_{max}$ as a function of the Reynolds number $Re$ and the Weber number $We$ (a) in linear scale and (b) in logarithmic scale. The parameters are $m=0.01$, $\eta =0.001$ and $n=5$.

The influence of the Weber number at fixed $Re=9000$ is isolated in figure 11 for three different azimuthal numbers. In all cases, $G_{max}$ is initially marginally below 1, followed by a sharp increase that leads to a local maximum which tends to have larger magnitude at smaller $n$. A further increase of the Weber number causes $G_{max}$ to decrease before remaining approximately constant for a comparatively large range of Weber numbers. Since the numerator of $G$ is proportional to $f$ and $\sigma$, and thus the inverse of the Weber number, this result implies that the magnitude of the surface distortion grows with $We$.

Figure 11. Energy gain $G_{max}$ as a function of the Weber number $We$ (a) in linear scale and (b) in logarithmic scale for $n=5$ (black), $n=6$ (blue) and $n=7$ (red). The critical Weber numbers $We_c$ corresponding to $G_{max}=1$ are marked by dashed lines. The parameters are $Re=9000$, $m=0.01$ and $\eta =0.001$.

Analogously, we investigate the scaling of $G_{max}$ with the Reynolds number at constant Weber number. The resulting $G_{max}$ for a density ratio $m=0.01$ and a viscosity ratio $\eta =0.001$ is presented in figures 12 and 13. We focus on two different regimes, the first relating to the local maximum of $G_{max}$ at relatively low Weber number seen in figure 11. The three considered cases are $n=5$ at $We=3000$, $n=6$ at $We=5000$ and $n=7$ at $We=6000$. The scaling of $G_{max}$ with the Reynolds number is presented in figure 12(a) for three azimuthal wavenumbers. Within the considered range of Reynolds numbers, $G_{max}$ linearly increases with $Re$. The slope of the lines nonetheless decreases with $n$. When pre-multiplied with the square of the azimuthal wavenumber $n$, the different lines effectively collapse (see figure 12b). The result thus establishes that the peak $G_{max}$ scales as

(5.1)\begin{equation} G_{max}\sim\frac{Re}{n^2}. \end{equation}

The second regime is related to the region of near-constant $G_{max}$ at relatively high $We$. For $We=20\,000$, $G_{max}$ again grows linearly with the Reynolds number for sufficiently large $Re$ (see figure 13a). As demonstrated in figure 13(b), the inverse proportionality of $G_{max}$ to the square of the azimuthal wavenumber does not, however, extend to the regime of higher Weber numbers.

Figure 12. Amplification ratio $G_{max}$ as a function of the Reynolds number $Re$ at the local maximum in terms of the Weber number: (a) $G_{max}$ and (b) $n^2 G_{max}$. Azimuthal wavenumbers $n=5$, $We \sim 3000$ (black); $n=6$, $We \sim 5000$ (blue); and $n=7$, $We \sim 6000$ (red). The parameters are $m=0.01$ and $\eta =0.001$.

Figure 13. Amplification ratio $G_{max}$ as a function of the Reynolds number $Re$ in the regime of relatively high Weber number: (a) $G_{max}$ and (b) $n^2 G_{max}$. Azimuthal wavenumbers $n=5$ (black), $n=6$ (blue) and $n=7$ (red). The parameters are $We=20\,000$, $m=0.01$ and $\eta =0.001$.

Overall, the results establish a linear scaling of $G_{max}$ with the Reynolds number throughout the considered range of Weber numbers. The differing scaling with the azimuthal wavenumber depending on the Weber number nonetheless points to differences in the evolution of the perturbation flow field. These differences are qualitatively illustrated by the perturbation fields in cross-planes at both the initial and final locations between the cases of low and high Weber number (see figures 14 and 15). At low Weber number, the highest perturbation amplitudes in the radial velocity component at the target location occur within the jet, at approximately two-thirds of the jet radius. In contrast, the most intense fluctuations at comparatively high $We$ are observed at the interface location.

Figure 14. Cross-planes showing (a) the optimal initial condition at $x_0=0$ and (b) its evolution at $x_1=11.8$. Vectors of the radial and azimuthal disturbance fields, $v'$ and $w'$, with colour contours of streamwise disturbance field, $u'$. The red line marks the position of the material interface. The parameters are $Re=7500$, $We=3000$, $m=0.01$, $\eta =0.001$, $n=5$ and $\omega =2.2$.

Figure 15. Cross-planes showing (a) the optimal initial condition at $x_0=0$ and (b) its evolution at $x_1=11.2$. Vectors of the radial and azimuthal disturbance fields, $v'$ and $w'$, with colour contours of streamwise disturbance field, $u'$. The red line marks the position of the material interface. The parameters are $Re=7500$, $We=10000$, $m=0.01$, $\eta =0.001$, $n=5$ and $\omega =1.7$.

We now turn to examining the effect of the viscosity ratio $m$ and the density ratio $\eta$ on the amplification mechanism. We maintain a tractable problem by keeping the remainder of the parameter space at constant values, $Re=6000$, $We=1000$ and $n=5$. Contours in figure 16 present the maximum surface tension energy gain $G_{max}$ for the range $0.01 \leqslant m \leqslant 0.335$ and $0.001 \leqslant \eta \leqslant 0.08$. We note that the viscosity ratio determines the ratio of the thicknesses of the shear layers in the two fluids, with higher values of $m$ corresponding to a relatively thicker shear layer in the gas phase. Moreover, the density ratio is related to the ratio of the perturbation kinetic energy in the two phases. The result shows that the growth of the interface distortion is highest for a combination of comparatively low viscosity ratio and high density ratio.

Figure 16. Contours of the energy gain $G_{max}$ as a function of the viscosity ratio $m$ and the density ratio $\eta$. The parameters are $Re=6000$, $We=1000$ and $n=5$.

The results presented in figures 10 and 11 established that for small Weber numbers, the objective functional remains marginally below unity before rising sharply once a certain threshold for $We$ is exceeded. In this context, we recall that the definition of the objective functional implies that for $G_{max}<1$, the perturbation surface tension energy at the target location is less than the total perturbation energy at the initial position, implying either passive decay or a trivial redistribution of initial perturbation kinetic energy into the interface. On the other hand, for cases with $G_{max}>1$, the perturbation interface energy at the target position exceeds the total perturbation energy at the initial position, which implies that energy has been extracted from the mean shear and converted into surface tension energy.

In the following, we seek to characterize the critical Weber number $We_c$ at which $G_{max}$ crosses unity, and which thus marks the boundary at which the distortion to the material interface grows in a non-trivial manner by means of a redistribution of energy from the mean shear. In terms of the energy budget, this implies that the right-hand side of (4.4) yields a positive value so that the contribution to the perturbation kinetic energy by the radial component at the target location $x_1$ has increased compared to the initial location $x_0$. The difference in the kinetic energy in the radial velocity component between the initial and target locations, $\varDelta \frac {1}{2}v'^2=\frac {1}{2}(v'^2(x_1)-v'^2(0))$, is presented in figure 18(a). The Weber numbers at which each of the curves crosses zero coincide with the critical Weber numbers identified in figure 11.

In the exploration of a potential scaling law for the critical Weber number, we turn to the critical Ohnesorge number,

(5.2)\begin{equation} Oh_c =\frac{\sqrt{We_c}}{Re}, \end{equation}

which relates the viscous forces to the inertial and surface tension forces at the critical Weber number. Figure 17 shows the logarithm of the Ohnesorge number, $\ln (Oh_c)$, for the corresponding critical Weber number as a function of the logarithm of the Reynolds number. A linear trend is observed with a slope of approximately $-3/2$, suggesting that $Oh_c\sim Re^{-3/2}$.

Figure 17. Scaling of the logarithm of the critical Ohnesorge number $Oh_c$ with the logarithm of the Reynolds number $Re$. Here $n=5$ (black), $n=6$ (blue) and $n=7$ (red). The parameters are $m=0.01$ and $\eta =0.001$. The black dotted line indicates a reference slope of $-3/2$.

Insight into the scaling of $Oh_c$ and thus $We_c$ is gained by extending the analysis of the energy budget presented in the preceding section, where it was shown that the PVC term is the major source of the surface tension energy. Since the measure for the gain of the interface energy $G$ relates the surface tension energy at a certain downstream location $x_1$ to the initial total energy of the perturbation, it follows that $G$ is proportional to the PVC term in (4.7). Previous results in figure 12(a) further showed that $G$ scales with $Re$. Thus,

(5.3)\begin{equation} \int_{\mathcal{\varOmega}} \frac{p'}{r}\frac{\partial }{\partial r}(rv')\,\text{d}\mathcal{\varOmega} \approx C_1 Re,\quad \text{with}\ C_1\geqslant0. \end{equation}

For the considered cases, $n^2\gg 1$, and also $n^2\gg \alpha ^2$. The surface tension energy term is thus approximated as

(5.4)\begin{equation} \int_{\mathcal{S}} \frac{1}{We} \left( - 1 + n^2 + \alpha^2 \right)f'v'\,\textrm{d}\mathcal{S} \approx \frac{C_2}{We}, \quad \text{with}\ C_2\geqslant0. \end{equation}

Using the energy budget in the radial velocity component (4.4), the condition at the critical Weber number of $\varDelta \frac {1}{2}v^2 \geqslant 0$ becomes

(5.5)\begin{align} &\int_{\varOmega}\frac{p'}{r}\frac{\partial }{\partial r}(rv')\,\textrm{d}\varOmega + {\int_{\mathcal{S}} \left( -\eta p'_g v'_g + \frac{2m}{Re}v'\frac{\partial v'_g}{\partial r} \right)\textrm{d}{\mathcal{S}}} \nonumber\\ &\quad - \int_{\mathcal{S}} {\frac{1}{Re}v'\frac{\partial v'}{\partial r}}\,\textrm{d}{\mathcal{S}} + \int_{\varOmega}{\frac{1}{Re} \frac{\partial^2}{\partial x \partial x} \left( \frac{1}{2}v'^2 \right)\textrm{d}\varOmega} + \int_{\mathcal{S}}{\frac{1}{We} \left( -1 + n^2 + \alpha^2 \right)f'v'}\,\textrm{d}{\mathcal{S}} \nonumber\\ &\quad -\int_{\varOmega} \frac{1}{Re} \left( \frac{\partial v'}{\partial r}\frac{\partial v'}{\partial r} + \frac{1}{r}\frac{\partial v'}{\partial \theta}\frac{1}{r}\frac{\partial v'}{\partial \theta} + \frac{\partial v'}{\partial x}\frac{\partial v'}{\partial x} + \frac{v'^2}{r^2} + \frac{2}{r^2}v'\frac{\partial w'}{\partial \theta} \right)\textrm{d}\varOmega \geqslant 0. \end{align}

Substitution of (5.3) and (5.4) into (5.5) yields

(5.6)\begin{equation} C_1 Re \geqslant \frac{C_2}{We} + \frac{C_3}{Re}, \end{equation}

with $C_3$ defined as

(5.7)\begin{align} C_3 &= \int_{\mathcal{S}} \left( \eta p'_g v'_g - 2m v'\frac{\partial v'_g}{\partial r} \right)\,\textrm{d}\mathcal{S} + \int_{\mathcal{S}} v'\frac{\partial v'}{\partial r}\,\textrm{d}\mathcal{S} - \int_{\varOmega} \frac{\partial^2}{\partial x \partial x} \left( \frac{1}{2}v'^2 \right)\,\textrm{d}\varOmega \nonumber\\ &\quad + \int_{\varOmega}\left( \frac{\partial v'}{\partial r}\frac{\partial v'}{\partial r} + \frac{1}{r}\frac{\partial v'}{\partial \theta}\frac{1}{r}\frac{\partial v'}{\partial \theta} + \frac{\partial v'}{\partial x}\frac{\partial v'}{\partial x} + \frac{v'^2}{r^2} + \frac{2}{r^2}v'\frac{\partial w'}{\partial \theta} \right) \textrm{d}\varOmega. \end{align}

The results in figure 8 showed that $C_3$ is small compared to the PVC and surface tension energy terms, especially in the early stage of the evolution, and thus $C_3 \ll C_1, C_2$. Hence, for $Re \gg 1$, we obtain

(5.8)\begin{equation} C_1 Re^3 \geqslant C_2 \frac{Re^2}{We_c}, \end{equation}

so that

(5.9)\begin{equation} Oh_c \geqslant C_4Re^{-{3}/{2}}, \end{equation}

where $C_4={C_2}/{C_1}$. The predicted exponent of $-3/2$ recovers the scaling observed in figure 17 and thus substantiates the power-law relation between the Reynolds number and the critical Ohnesorge number. Substitution of (5.2) yields the associated critical Weber number:

(5.10)\begin{equation} We_c \geqslant \frac{C_4^2}{Re}. \end{equation}

For Weber numbers below the critical value, $We<We_c$, the energy gain attains values less than unity, and the surface tension energy gain may thus be attributed to a redistribution of the initial perturbation kinetic energy without meaningful extraction of energy from the mean shear. In this regime, the target location $x_1$ at which $G_{max}$ is attained is expected to be comparatively close to the initial position $x=0$ so as to minimize any losses of perturbation kinetic energy due to viscous dissipation. Results presented in figure 18(b) showing the target location $x_1$ as a function of Weber number confirm this hypothesis. For all considered azimuthal wavenumbers, the critical Weber numbers $We_c$ coincide with an abrupt downstream shift of the target locations.

Figure 18. (a) Perturbation kinetic energy gain in the radial direction, $\varDelta \frac {1}{2}v^2$, as a function of Weber number $We$. (b) Target location $x_1$ as a function of Weber number. Azimuthal wavenumbers $n=5$ (black), $n=6$ (blue) and $n=7$ (red). The parameters are $Re=9000$, $m=0.01$ and $\eta =0.001$. Dashed lines mark the critical Weber numbers $We_c$.

6. Nonlinear simulation

So far, we have demonstrated the algebraic mechanism for generating interface distortion exclusively by means of linear analyses. In this setting, the nonlinear higher-order terms in the governing equation are neglected which in turn allows the identification of optimal solutions in the interface energy norm as combinations of the non-orthogonal eigenfunctions. In the following, we demonstrate the realizability of the identified mechanism in a nonlinear setting by means of a simulation of the full two-fluid Navier–Stokes equations.

The nonlinear simulations employ the Basilisk solver (Popinet Reference Popinet2015, Reference Popinet2018) which is based on a conservative, non-diffusive geometric volume-of-fluid scheme to track the material interface. The evolution of the interface is captured by advecting a volume-fraction function coupled with the incompressible Navier–Stokes equations. An orthogonally structured grid of size $[N_x, N_y, N_z]=[1536\times 256\times 256]$ and spatial extent $[L_x, L_y, L_z]=[60R, 10R, 10R]$ is used. In this setting, the radius of the jet is defined as $r = (y^2+z^2)^{1/2}$ in the $y$–$z$ plane. A homogeneous uniform grid spacing of $\Delta x=0.004$ is applied in all dimensions. We note that the comparatively large streamwise extent of the computational domain prevents any interaction of the jet with the outflow plane during the runtime of the simulations.

In a first step, a base flow is computed by imposing the profile (2.2) at the inflow of the computational domain. After the initial transient has cleared, the base flow profile is superimposed with the optimal initial condition identified in the solution of the optimization problem (3.1) for the parameters $Re=7784$, $We=4865$, $m=0.036$, $\eta =0.0012$, $n=6$ and $\omega =2.5$. Homogeneous Neumann boundary conditions are imposed on both the side walls and the outlet:

(6.1)\begin{align} \left.\begin{array}{ll@{}} \dfrac{\partial\boldsymbol{u}\boldsymbol{\cdot}\hat{\boldsymbol{e}}_n}{\partial n} = 0 & {\rm if}\ \dfrac{\partial\boldsymbol{u}\boldsymbol{\cdot}\hat{\boldsymbol{e}}_n}{\partial n} > 0,\\ \boldsymbol{u}\boldsymbol{\cdot}\hat{\boldsymbol{e}}_n = 0 & \text{otherwise}, \end{array}\right\}\end{align}

where $\hat {\boldsymbol {e}}_n$ is the face-normal vector pointing outward of the domain. The convergence of the simulation results was verified by increasing the resolution as well as changing the domain size in the $y$ and $z$ directions.

An isosurface of the material interface at time $t=0.24$ is presented in figure 19(a). Colours indicate the local amplitude of the interface distortion. The linear result, rescaled to the same initial amplitude of the perturbations, is presented in figure 19(b) and shows excellent qualitative agreement.

Figure 19. Isosurfaces of the material interface, coloured by the local displacement of the interface. (a) Linear stability theory. (b) Nonlinear simulation based on a geometric volume-of-fluid formulation. The parameters are $Re=7784$, $We=4865$, $m=0.036$, $\eta =0.0012$, $n=6$ and $\omega =2.5$.

The envelopes of the amplitudes of the interface perturbations computed in linear analysis and the nonlinear simulations are plotted in figure 20. The nonlinear simulation qualitatively recovers the non-monotonic growth of the perturbation observed in the linear analysis. In both cases, the global maximum is preceded by two local maxima past which the amplitude of the interface distortion decays before amplifying again. The prediction by linear theory and the nonlinear simulation are in excellent agreement in the vicinity of the nozzle. With growing distance to the exit of the nozzle, a moderate discrepancy between the two results appears. The peak amplitude of the interface distortion is recorded at $x = 11.2$ in the linear analysis while it is located at $x = 12.37$ in the nonlinear simulation. The magnitude of the maximum interface distortion also attains a moderately lower value of $|\,f(x)/\,f_0| = 16.05$ in the nonlinear simulation, compared to a value of $|\,f(x)/\,f_0| = 16.81$ predicted in the linear analysis. This disparity can be attributed to the relaxation of the base flow in the nonlinear simulation which leads to a thickening of the shear layer with increasing distance to the inflow location at $x=0$. In contrast, the linear analysis assumes a constant base flow profile which matches that at the inflow location of the nonlinear simulation. An evaluation of the surface tension energy gain $G_{max}$ for the relaxed base flow at $x = 10$ results in a reduced amplification ratio, $|\,f(x)/\,f_0| = 14.08$. Overall, the qualitative agreement between linear analysis and nonlinear simulations clearly substantiates the realizability of the presented mechanism.

Figure 20. Envelope of the interface distortion of the liquid jet. Linear stability theory (black) and nonlinear simulation based on a geometric volume-of-fluid formulation (blue). The parameters are $Re=7784$, $We=4865$, $m=0.036$, $\eta =0.0012$, $n=6$ and $\omega =2.5$.