1 Introduction
Coaxial jets are often employed as an effective way to rapidly mix two different fluid streams, such as in industrial burners and airblast atomisers, to mention just two examples. In addition to their mixing properties, these flows are also well known for the possibility of reducing noise generation compared to the single jet configuration (Williams, Ali & Anderson Reference Williams, Ali and Anderson1969), resulting in several aerospace applications. Both noise and mixing properties strongly depend on the dynamics of the large-scale vortical structures produced by the inherent instabilities of the flow field, whose deep understanding is therefore relevant not only from the academic viewpoint.
The incompressible flow produced by two coaxial jets is characterised by several non-dimensional parameters that are associated with the nozzle geometry and the characteristics of the two incoming fluid streams: for instance the outer-to-inner pipe diameter ratio,
$R_{D}$
, the ratio between the maximum velocity of the outer and inner jets,
$R_{u}$
, the boundary layer thickness at the nozzle exit and the free stream turbulence. Past experimental and numerical investigations (see for example Ko & Kwan Reference Ko and Kwan1976; Dahm, Clifford & Tryggvanson Reference Dahm, Clifford and Tryggvanson1992; Rehab, Villermaux & Hopfinger Reference Rehab, Villermaux and Hopfinger1997; Villermaux & Rehab Reference Villermaux and Rehab2000; da Silva, Balarac & Metais Reference da Silva, Balarac and Metais2003; Segalini & Talamelli Reference Segalini and Talamelli2011) have established the central role of
$R_{u}$
in the coupling between the inner and outer shear layer instabilities. For
$R_{u}\gg 1$
a strong interaction occurs which results in the synchronised roll up of both the shear layers, the so-called lock-in phenomenon. In this regime the coaxial jet dynamics is driven by the Kelvin–Helmholtz instability of the outer shear layer (da Silva et al.
Reference da Silva, Balarac and Metais2003). Conversely for
$R_{u}\ll 1$
, a weak interaction is observed where the outer stream mainly acts as a ‘coflow’ (Djeridane Reference Djeridane1994), without substantially modifying either the inner jet dynamics or the main flow instability, which is now driven by the inner jet shear layer (Segalini Reference Segalini2010). Within the range of nearly unitary velocity ratios, the development of vortical structures in the flow is strongly affected by the presence or absence of a thick duct wall separating the two streams (Buresti, Petagna & Talamelli Reference Buresti, Petagna and Talamelli1994; Segalini Reference Segalini2010). Indeed, in the former case, the near flow field is dominated by an alternate vortex shedding from the blunt edge of the separating wall, which considerably enhances the mixing and entrainment processes of the two jets (Talamelli et al.
Reference Talamelli, Segalini, Orlu and Buresti2013). The existence of three different instability regimes has also been described in the experimental work by Segalini & Talamelli (Reference Segalini and Talamelli2011), based on an extensive analysis of the dominant frequencies characterising the unsteady flow field. The authors have also provided an effective relationship for the proper scaling of the shedding frequency in the intermediate range of velocity ratios,
$0.75\lesssim R_{u}\lesssim 1.6$
, improving the previous one proposed by Buresti et al. (Reference Buresti, Petagna and Talamelli1994).
Most of the studies addressing the coaxial jet stability properties have been performed based on a local approach, using analytical and experimental velocity profiles (Michalke Reference Michalke1993; Juniper & Candel Reference Juniper and Candel2003; Talamelli & Gavarini Reference Talamelli and Gavarini2006; Gloor, Obrist & Kleiser Reference Gloor, Obrist and Kleiser2013). Within this framework, Talamelli & Gavarini (Reference Talamelli and Gavarini2006) have investigated the inviscid stability properties of a family of axial velocity profiles, varying several related parameters. In particular, the authors have shown that an absolute axisymmetric instability occurs, and that it depends on the local minimum velocity
$\widetilde{U}_{w}$
in the backflow region representing the duct wall wake. The existence of this instability is also restricted to specific values of the velocity ratio, and it was found that the range of
$R_{u}$
where the instability occurs shrinks around
$R_{u}\approx 1.3$
as
$\widetilde{U}_{w}$
is increased, i.e. while moving downstream of the nozzle. Both the lower
$R_{u}$
instability threshold and the absolute instability frequency, approximately equal to 0.5, are found to weakly depend on the momentum thickness of the boundary layer. A quasi-linear dependence of the mode frequency on
$R_{u}$
is reported in the same work, and a good agreement between numerical and experimental values of the non-dimensional frequency is obtained when introducing a scaling based on the average velocity of the two streams and on the duct wall wake thickness.
Besides the inherent convective nature of the jet flow instabilities, the existence of a pocket of absolute instability and the annular vortex shedding observed in the experiments suggest that, for certain values of
$Re$
and
$R_{u}$
, the instability of two coaxial jets separated by a thick wall can be ascribed to the onset of an unstable global mode. This mechanism, which leads to a Hopf bifurcation of the steady axisymmetric base flow, is expected to be dominant in the range of unitary velocity ratios, implying a ‘transition’ from a noise-amplifier to a fluid-oscillator behaviour (Huerre & Rossi Reference Huerre, Rossi, Godrèche and Manneville1998). As pointed out by Talamelli et al. (Reference Talamelli, Segalini, Orlu and Buresti2013), the onset of self-excited oscillations of the whole flow field is of great interest, since it provides a ‘natural’ forcing mechanism of the flow, increasing the mixing without requiring the introduction of additional energy in the fluid system, such as in the case of an active control device. At the same time, this self-excited instability can become a source of loud noise (Olsen & Karchmer Reference Olsen and Karchmer1976). These observations motivate us to investigate the stability properties of the coaxial jet flow using a global approach, thus accounting for strong non-parallel effects in the near flow field caused by the thick separating wall. Both modal and non-modal analyses are performed to characterise the different instabilities of the flow in the range of
$0.5\lesssim R_{u}\lesssim 2$
and
$1000\lesssim Re\lesssim 5000$
. On the one hand, the non-modal analysis is more appropriate to describe the relevant instability mechanisms when a noise-amplifier flow behaviour occurs (Garnaud et al.
Reference Garnaud, Lesshafft, Schmid and Huerre2013). On the other hand, the modal analysis is extended to consider the structural sensitivity of the expected unstable global mode with respect to linearised perturbations. This analysis is performed to identify those flow regions which are responsible for the development of self-excited flow oscillations, as done by Giannetti & Luchini (Reference Giannetti and Luchini2007) for the cylinder wake. Finally, both axisymmetric and fully three-dimensional direct numerical simulations (DNS) of the flow are performed to assess the validity and relevance of the stability results.
The paper is organised as follows. The flow configuration and the governing equations are introduced in § 2 where the different approaches employed to characterise the flow stability properties are shortly recalled. Some details about their numerical implementation are given in § 3, along with validation tests on the single round jet configuration. The steady axisymmetric base flow is described in § 4. Results obtained from the modal stability analysis are discussed in § 5. More precisely, § 5.1 presents the global eigenspectra computed for different perturbation wavenumbers, with the identification of the unstable global mode responsible for the onset of the vortex shedding in the range of nearly unitary velocity ratios. The corresponding neutral curve is tracked in the parameter plane
$Re$
–
$R_{u}$
, thus defining the domain of linear asymptotic instability of the axisymmetric base flow. The leading direct and adjoint eigenfunctions are illustrated in § 5.2 together with the associated wavemaker region. Then the transient dynamics of the linearised flow perturbations is investigated in § 5.3. Finally, the results obtained from the direct numerical simulations of the flow are compared with the modal stability results in § 6 and some conclusions are drawn in § 7.
2 Mathematical formulation
2.1 Flow configuration
Figure 1. Sketch of the computational domain
$\unicode[STIX]{x1D6FA}_{c}$
(azimuthal plane) employed for the investigation of the coaxial jet flow. All lengths are made non-dimensional using the inner pipe diameter
$\widetilde{D}_{i}$
. The black thick line is used to represent the solid wall boundary
$\unicode[STIX]{x1D6E4}_{w}$
while the blue shaded area denotes the region of mesh refinement. The outer-to-inner pipe diameter ratio is fixed to
$R_{D}=2$
, while the thickness of the duct wall is 0.1.
In the present paper we are concerned with the incompressible flow produced by two coaxial jets in an unbounded region of quiescent fluid. The axisymmetric computational domain
$\unicode[STIX]{x1D6FA}_{c}$
and the adopted cylindrical coordinate system
$(r,\unicode[STIX]{x1D703},z)$
are illustrated in figure 1. All lengths are made non-dimensional using the inner pipe diameter
$\widetilde{D}_{i}$
. The flow configuration is chosen to be very similar to that investigated by Segalini & Talamelli (Reference Segalini and Talamelli2011), with an outer-to-inner pipe diameter ratio of
$R_{D}=2$
and a non-dimensional thickness of the duct wall of
$s=0.1$
. The main difference with respect to the geometry employed by these authors consists in the presence of the additional solid wall on the plane
$z=0$
, which surrounds the annular pipe. The flow dynamics is described by the unsteady incompressible Navier–Stokes equations which are made dimensionless using
$\widetilde{D}_{i}$
, the maximum velocity
$\widetilde{U}_{i}$
of the inner fluid stream at the pipe inlet and the constant density
$\tilde{\unicode[STIX]{x1D70C}}$
:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{U}}{\unicode[STIX]{x2202}t}+(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{U}+\unicode[STIX]{x1D735}P-\displaystyle \frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{U}=\mathbf{0},\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}=0.\end{array}\right\}\end{eqnarray}$$
The Reynolds number
$Re$
is defined as
$Re=\widetilde{U}_{i}\widetilde{D}_{i}/\tilde{\unicode[STIX]{x1D708}}$
,
$\tilde{\unicode[STIX]{x1D708}}$
being the kinematic viscosity of the fluid. Therefore, the total flow state
$\boldsymbol{Q}=\{\boldsymbol{U},P\}^{\text{T}}$
is specified by the velocity field
$\boldsymbol{U}=\boldsymbol{U}(r,\unicode[STIX]{x1D703},z,t)$
with components
$(U_{r},U_{\unicode[STIX]{x1D703}},U_{z})^{\text{T}}$
, and by the reduced pressure field
$P=P(r,\unicode[STIX]{x1D703},z,t)$
. With reference to figure 1, no-slip conditions are imposed on the wall boundary
$\unicode[STIX]{x1D6E4}_{w}$
while the following condition is assumed on both
$\unicode[STIX]{x1D6E4}_{t}$
and
$\unicode[STIX]{x1D6E4}_{o}$
:
$$\begin{eqnarray}\displaystyle \frac{1}{Re}\frac{\unicode[STIX]{x2202}\boldsymbol{U}}{\unicode[STIX]{x2202}\hat{\boldsymbol{n}}}-P\hat{\boldsymbol{n}}=\mathbf{0},\end{eqnarray}$$
where
$\hat{\boldsymbol{n}}$
denotes the outward unit normal vector. On the inlet boundaries
$\unicode[STIX]{x1D6E4}_{in,i}$
and
$\unicode[STIX]{x1D6E4}_{in,o}$
, the following Dirichlet conditions hold:
$$\begin{eqnarray}\boldsymbol{U}|_{\unicode[STIX]{x1D6E4}_{in,i}}=U_{z,i}(r)\hat{\boldsymbol{z}},\quad \boldsymbol{U}|_{\unicode[STIX]{x1D6E4}_{in,o}}=U_{z,o}(r)\hat{\boldsymbol{z}},\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle & \displaystyle U_{z,i}(r)=\tanh \left[b_{i}\left(1-2r\right)\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle U_{z,o}(r)=R_{u}\tanh \left[b_{o}\left(1-\left|\displaystyle \frac{2r-(R_{o,1}+R_{o,2})}{R_{o,1}-R_{o,2}}\right|\right)\right], & \displaystyle\end{eqnarray}$$
$R_{o,1}$
and
$R_{o,2}$
denote the internal and external non-dimensional radii of the annular duct and
$R_{u}$
represents the velocity ratio between the two jets, defined as
$R_{u}=\max U_{z,o}/\text{max}\,U_{z,i}$
. The two above analytical expressions are chosen to reproduce the experimental velocity profiles reported in the work by Segalini (Reference Segalini2010). More precisely, the two parameters
$b_{i}$
and
$b_{o}$
are related to the thickness of the boundary layers within the nozzle. In the present analysis both values are fixed to
$5$
. Indeed past works have shown that the pipe boundary layer features have only a weak influence on the flow stability properties, especially in the vortex shedding regime (Talamelli & Gavarini Reference Talamelli and Gavarini2006). The incompressible Navier–Stokes equations with the aforementioned boundary conditions are solved numerically by Newton’s method to compute the base flow
$\boldsymbol{Q}_{b}=\{\boldsymbol{U}_{b}(r,z),P_{b}(r,z)\}^{\text{T}}$
on top of which the linear stability analysis is performed, as detailed in the next paragraph. In this respect, the base flow represents the numerical approximation of an exact solution of the Navier–Stokes equations, as done for the single jet by Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013).
2.2 Linear global stability
As already mentioned, the onset of the various flow instabilities is investigated by performing a modal and a non-modal analysis of the linearised perturbations evolving on top of the steady axisymmetric solution to (2.1),
$\boldsymbol{Q}_{b}$
. The axisymmetric hypothesis is no longer valid for the total linearised perturbation field
$\boldsymbol{q}=\{\boldsymbol{u},p\}^{\text{T}}$
which is expanded in Fourier modes of azimuthal wavenumber
$m\in \mathbb{N}$
,
$$\begin{eqnarray}\boldsymbol{q}(r,\unicode[STIX]{x1D703},z,t)=\mathop{\sum }_{m=-\infty }^{\infty }\boldsymbol{q}_{m}(r,z,t)\text{e}^{\text{i}m\unicode[STIX]{x1D703}}.\end{eqnarray}$$
The complex-conjugate symmetry
$\boldsymbol{q}_{m}=\boldsymbol{q}_{-m}^{\ast }$
holds for these modes, and the notation
$(\cdot )^{\ast }$
is employed here and in the following to indicate the complex-conjugate of a given complex-valued quantity. At each wavenumber
$m$
, the corresponding perturbation field
$\boldsymbol{q}_{m}=\{\boldsymbol{u}_{m},p_{m}\}^{\text{T}}$
is governed by the linear dynamical system
$$\begin{eqnarray}{\mathcal{B}}\frac{\unicode[STIX]{x2202}\boldsymbol{q}_{m}}{\unicode[STIX]{x2202}t}={\mathcal{L}}_{m}\boldsymbol{q}_{m},\end{eqnarray}$$
where the operators
${\mathcal{B}}$
and
${\mathcal{L}}_{m}$
are expressed as follows:
$$\begin{eqnarray}{\mathcal{L}}_{m}\boldsymbol{q}_{m}\equiv \left(\begin{array}{@{}c@{}}-{\mathcal{C}}_{m}(\boldsymbol{u}_{m},\boldsymbol{U}_{b})+Re^{-1}\unicode[STIX]{x1D6FB}_{m}^{2}\boldsymbol{u}_{m}-\unicode[STIX]{x1D735}_{m}p_{m}\\ \unicode[STIX]{x1D735}_{m}\boldsymbol{\cdot }\boldsymbol{u}_{m}\end{array}\right),\quad {\mathcal{B}}\frac{\unicode[STIX]{x2202}\boldsymbol{q}_{m}}{\unicode[STIX]{x2202}t}\equiv \left(\begin{array}{@{}c@{}}\unicode[STIX]{x2202}\boldsymbol{u}_{m}/\unicode[STIX]{x2202}t\\ 0\end{array}\right),\end{eqnarray}$$
with
$$\begin{eqnarray}{\mathcal{C}}_{m}(\boldsymbol{u}_{m},\boldsymbol{U}_{b})=(\boldsymbol{U}_{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{m})\boldsymbol{u}_{m}+(\boldsymbol{u}_{m}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{0})\boldsymbol{U}_{b}.\end{eqnarray}$$
The notation
$\unicode[STIX]{x1D735}_{m}(\cdot )$
,
$\unicode[STIX]{x1D735}_{m}\boldsymbol{\cdot }(\cdot )$
and
$\unicode[STIX]{x1D6FB}_{m}^{2}(\cdot )$
is used to indicate the Fourier transformed gradient, divergence and vector Laplacian operators, respectively. For the linearised flow field the same boundary conditions introduced in § 2.1 are applied with homogeneous data, as follows from the linearisation of the original nonlinear differential problem.
2.2.1 Global modes and wavemaker analysis
For each wavenumber
$m$
the modal stability analysis of the axisymmetric base flow is performed by casting the perturbation field
$\boldsymbol{q}_{m}$
in the normal mode form
$\boldsymbol{q}_{m}(r,z,t)=\hat{\boldsymbol{q}}_{m}(r,z)\text{e}^{\unicode[STIX]{x1D706}t}$
. This ansatz leads to the following generalised eigenvalue problem for the global mode
$\hat{\boldsymbol{q}}_{m}(r,z)$
and the associated eigenvalue
$\unicode[STIX]{x1D706}=\unicode[STIX]{x1D70E}+\text{i}\unicode[STIX]{x1D714}\in \mathbb{C}$
:
$$\begin{eqnarray}\unicode[STIX]{x1D706}{\mathcal{B}}\hat{\boldsymbol{q}}_{m}={\mathcal{L}}_{m}\hat{\boldsymbol{q}}_{m},\end{eqnarray}$$
with
$\unicode[STIX]{x1D70E}$
and
$\unicode[STIX]{x1D714}$
denoting the real and imaginary part of
$\unicode[STIX]{x1D706}$
, respectively. A global instability clearly arises when
$\unicode[STIX]{x1D70E}>0$
for some values of the governing parameters
$Re$
and
$R_{u}$
. The corresponding adjoint global modes
$\boldsymbol{q}_{m}^{\dagger }(r,z,t)=\hat{\boldsymbol{q}}_{m}^{\dagger }(r,z)\text{e}^{-\unicode[STIX]{x1D706}^{\ast }t}$
are introduced as well with
$$\begin{eqnarray}\unicode[STIX]{x1D706}^{\ast }{\mathcal{B}}\hat{\boldsymbol{q}}_{m}^{\dagger }={\mathcal{L}}_{m}^{\dagger }\hat{\boldsymbol{q}}_{m}^{\dagger },\end{eqnarray}$$
where
${\mathcal{L}}_{m}^{\dagger }$
is the adjoint operator of
${\mathcal{L}}_{m}$
with respect to the Hermitian scalar product
$$\begin{eqnarray}\langle \boldsymbol{q}_{m}^{\dagger },\boldsymbol{q}_{m}\rangle _{\unicode[STIX]{x1D6FA}_{c}}=\int _{\unicode[STIX]{x1D6FA}_{c}}\left(\boldsymbol{q}_{m}^{\dagger }\right)^{\ast }\boldsymbol{\cdot }\boldsymbol{q}_{m}r\,\text{d}r\,\text{d}z.\end{eqnarray}$$
In particular, when a global mode becomes unstable, leading to an oscillator-type behaviour, the knowledge of both the direct and the adjoint eigenfunctions allows one to identify the ‘core’ region of the self-excited instability mechanism (Giannetti & Luchini Reference Giannetti and Luchini2007). This is achieved by performing a structural sensitivity analysis of the eigenvalue problem (2.9) with respect to a ‘localised’ structural perturbation of the operator
${\mathcal{L}}_{m}$
in the form of a localised feedback from the velocity field:
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}{\mathcal{L}}_{m}=\displaystyle \frac{1}{r}\unicode[STIX]{x1D6FF}(r-r_{0},z-z_{0})\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D63E}_{0} & 0\\ 0 & 0\end{array}\right).\end{eqnarray}$$
Here
$\unicode[STIX]{x1D63E}_{0}$
is a constant feedback tensor,
$(r_{0},z_{0})$
denote the coordinates of the point in the azimuthal plane where the feedback acts and
$\unicode[STIX]{x1D6FF}(r-r_{0},z-z_{0})$
is the two-dimensional Dirac delta function. It is worthwhile to note that, consistently with the axisymmetric nature of the underlying base flow, the considered structural perturbation is localised in
$r$
and
$z$
but not in
$\unicode[STIX]{x1D703}$
, with the wavemaker region resulting axisymmetric. By employing the formulation adopted by Pralits, Brandt & Giannetti (Reference Pralits, Brandt and Giannetti2010), the structural sensitivity analysis leads eventually to the definition of a sensitivity tensor field
$\unicode[STIX]{x1D64E}_{m}(r,z)$
for each considered wavenumber
$m$
$$\begin{eqnarray}\unicode[STIX]{x1D64E}_{m}(r,z)=\displaystyle \frac{(\hat{\boldsymbol{u}}_{m}^{\dagger })^{\ast }\otimes \hat{\boldsymbol{u}}_{m}}{\langle \hat{\boldsymbol{q}}_{m}^{\dagger },{\mathcal{B}}\hat{\boldsymbol{q}}_{m}\rangle _{\unicode[STIX]{x1D6FA}_{c}}},\end{eqnarray}$$
where
$\otimes$
denotes the dyadic product between two vector fields. A corresponding scalar sensitivity map in the azimuthal plane
$(r,z)$
is easily obtained by plotting at each spatial point a suitable norm of the tensor, such as, for instance, the Frobenius norm,
$\Vert \unicode[STIX]{x1D64E}_{m}(r,z)\Vert _{F}$
, which will be adopted in the present study. Other norm definitions can be considered as well, since the wavemaker identification has been shown not to depend on this particular choice (Camarri & Giannetti Reference Camarri and Giannetti2010; Carini, Giannetti & Auteri Reference Carini, Giannetti and Auteri2014b
).
2.2.2 Transient growth
In contrast with an oscillator-type dynamics, when the flow is dominated by a noise-amplifier behaviour, a non-modal approach allows the characterisation of those instabilities which are promoted by a relevant transient amplification of the linearised perturbations, due to non-normal interactions among the stable eigenmodes. This problem is commonly addressed by looking at the most amplified initial condition
$\boldsymbol{u}_{m,0}=\boldsymbol{u}_{m}(r,z,0)$
over a finite time horizon
$\unicode[STIX]{x1D70F}$
, i.e. the optimal perturbation for
$\unicode[STIX]{x1D70F}$
(Reddy & Henningson Reference Reddy and Henningson1993; Schmid Reference Schmid2007). The optimal amplification factor for the perturbation wavenumber
$m$
and a given time horizon
$\unicode[STIX]{x1D70F}$
,
$G_{m}^{opt}(\unicode[STIX]{x1D70F})$
, is defined as
$$\begin{eqnarray}G_{m}^{opt}(\unicode[STIX]{x1D70F})=\max _{\boldsymbol{ u}_{m,0}}\displaystyle \frac{\Vert \boldsymbol{u}_{m}(r,z,\unicode[STIX]{x1D70F})\Vert }{\Vert \boldsymbol{u}_{m,0}\Vert },\end{eqnarray}$$
with
$\Vert \boldsymbol{u}_{m}\Vert ^{2}=\int _{\unicode[STIX]{x1D6FA}_{c}}|\boldsymbol{u}_{m}|^{2}r\,\text{d}r\,\text{d}z$
. The velocity field at
$t=\unicode[STIX]{x1D70F}$
can be formally expressed in terms of the initial condition
$\boldsymbol{u}_{m,0}$
by introducing the propagator
${\mathcal{P}}_{m}(\unicode[STIX]{x1D70F})$
associated with the linearised system (2.6)
$$\begin{eqnarray}\boldsymbol{u}_{m}|_{t=\unicode[STIX]{x1D70F}}={\mathcal{B}}_{r}{\mathcal{P}}_{m}(\unicode[STIX]{x1D70F})\boldsymbol{q}_{m,0}={\mathcal{B}}_{r}\text{e}^{{\mathcal{L}}_{m}\unicode[STIX]{x1D70F}}\boldsymbol{q}_{m,0}={\mathcal{B}}_{r}\text{e}^{{\mathcal{L}}_{m}\unicode[STIX]{x1D70F}}{\mathcal{B}}_{p}\boldsymbol{u}_{m,0},\end{eqnarray}$$
where
${\mathcal{B}}_{r}$
and
${\mathcal{B}}_{p}$
denote the restriction and prolongation operators between the total flow state
$\boldsymbol{q}_{m}$
and the associated velocity field ‘component’
$\boldsymbol{u}_{m}$
:
$$\begin{eqnarray}{\mathcal{B}}_{r}\boldsymbol{q}_{m}=\boldsymbol{u}_{m},\quad {\mathcal{B}}_{p}\boldsymbol{u}_{m}=\left(\begin{array}{@{}c@{}}\boldsymbol{u}_{m}\\ 0\end{array}\right).\end{eqnarray}$$
Both
$G_{m}^{opt}(\unicode[STIX]{x1D70F})$
and the optimal perturbation
$\boldsymbol{u}_{m,0}^{opt}$
can be obtained as the leading eigenpair of the following symmetric positive-definite eigenvalue problem
$$\begin{eqnarray}{\mathcal{B}}_{r}{\mathcal{P}}_{m}^{\dagger }(\unicode[STIX]{x1D70F}){\mathcal{P}}_{m}(\unicode[STIX]{x1D70F}){\mathcal{B}}_{p}\boldsymbol{u}_{m,0}^{opt}=G_{m}^{opt}(\unicode[STIX]{x1D70F})^{2}\boldsymbol{u}_{m,0}^{opt},\end{eqnarray}$$
where
${\mathcal{P}}_{m}^{\dagger }(\unicode[STIX]{x1D70F})$
is the adjoint propagator of
${\mathcal{P}}_{m}(\unicode[STIX]{x1D70F})$
, i.e.
${\mathcal{P}}_{m}^{\dagger }(\unicode[STIX]{x1D70F})=\text{e}^{-{\mathcal{L}}_{m}^{\dagger }\unicode[STIX]{x1D70F}}$
.
3 Numerical methods
The steady solutions and the stability problems are numerically solved with PaStA, a Fortran90 code written in primitive variables and based on the finite element method (Canton Reference Canton2013; Canton, Schlatter & Örlü Reference Canton, Schlatter and Örlü2016). The equations are discretised on an unstructured, two-dimensional mesh composed of triangular elements, corresponding to an azimuthal plane of the considered cylindrical domain. Standard Taylor–Hood
$\mathbb{P}_{2}$
–
$\mathbb{P}_{1}$
finite elements are employed for the unknown velocity and pressure fields. In our numerical procedure the axisymmetric base flow is first computed by means of Newton iterations and then, for each considered wavenumber
$m$
, the corresponding global modes and optimal perturbations are extracted making use of the implicitly restarted Arnoldi algorithm implemented in the ARPACK library (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1998). More precisely a shift-invert spectral transformation is applied to the eigenproblem (2.9), with the Jacobian matrix being fully assembled. In contrast, a time stepper approach (Barkley, Blackburn & Sherwin Reference Barkley, Blackburn and Sherwin2008) is adopted for the computation of the optimal perturbations, through direct-adjoint iterations coupled to the aforementioned Arnoldi algorithm. The Crank–Nicolson time-scheme is used for the time integration of both the direct and the adjoint systems. In particular, a discrete adjoint approach at the spatial level is adopted, and the proper boundary conditions for the adjoint problem are thus accounted for automatically. The LOCA library of continuation algorithms (Salinger et al.
Reference Salinger, Bou-Rabee, Pawlowski, Wilkes, Burroughs, Lehoucq and Romero2002) is used to track the neutral curve associated with the leading eigenmode in the parameter plane
$Re$
–
$R_{u}$
. All matrix inversions are handled using the sparse LU solvers provided with the software library MUMPS (Amestoy, Duff & L’Excellent Reference Amestoy, Duff and L’Excellent2000).
In order to assess the accuracy of our stability computations with respect to the adopted discretisation, eleven meshes, characterised by different spatial extents, have been employed. The main features of these meshes are reported in table 1. For all the grids, a relevant portion of the inlet pipes has been modelled, with
$z_{i}\geqslant 5$
, in order to ensure the correct approximation of the physical adjoint eigenfunctions and of the optimal perturbations, as recommended by Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013). In the rest of the paper, we will present and discuss the results obtained using grids
$M_{60}$
–
$M_{100,z_{i}15}$
. A sensitivity analysis of these results with respect to the domain extents is reported in appendix A. All the meshes are characterised by a grid size
$\unicode[STIX]{x0394}\ell$
varying from
$\unicode[STIX]{x0394}\ell =0.4$
close to
$\unicode[STIX]{x1D6E4}_{o}$
, down to a minimum value of
$\unicode[STIX]{x0394}\ell =0.009$
within the refinement zone (blue shaded area in figure 1) these sizes of the elements have been determined following a mesh refinement independence study.
Finally, nonlinear DNS of the considered flow have been carried out by means of the spectral element code nek5000 (Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2008). As in the case of the finite element method, a classical Galerkin approximation is used to spatially discretise the governing equations. In this case, the velocity and the pressure spaces are spanned by Lagrange polynomial interpolants based on tensor-product arrays of Gauss–Lobatto–Legendre quadrature points in each local element. The velocity polynomial degree is
$N$
, two times higher than that employed for the pressure (
$\mathbb{P}_{N}$
–
$\mathbb{P}_{N-2}$
formulation); in the present computations
$N=6$
is chosen as a good compromise between accuracy and computational cost. Time integration employs an explicit second-order extrapolation for the advection terms and an implicit second-order backward differentiation for the viscous term. Both axisymmetric and fully three-dimensional simulations are performed. For the axisymmetric case the computational grid consists of
$3920$
elements, while for the three-dimensional case the total number of elements is equal to
$64\,160$
with
$32$
planes being employed in the azimuthal direction.
Table 1. Characteristics of the meshes used in the present study. All the meshes are named
$M_{z_{max}}$
so as to quickly recall the most important mesh parameter; where duplicates are present the changing parameter is reported as a second subscript.
$N_{e}$
denotes the total number of triangular elements while
$N_{dofs}$
denotes the corresponding total number of degrees of freedom. All the meshes are characterised by a grid size
$\unicode[STIX]{x0394}\ell$
varying from
$\unicode[STIX]{x0394}\ell _{max}=0.4$
close to the outflow border
$\unicode[STIX]{x1D6E4}_{o}$
, down to a minimum value of
$\unicode[STIX]{x0394}\ell _{min}=0.009$
within the refinement zone (blue shaded area in figure 1); these sizes of the elements have been determined following a mesh refinement independence study. See figure 1 for the definition of the other parameters.
3.1 Validation case: the single jet
Figure 2. Sketch of the computational domain employed for the stability analysis of the single jet configuration.
The global stability properties of the flow produced by an incompressible, round viscous jet have been carefully described by Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013). The same flow configuration is considered here as a test case for our numerical set-up. Similarly to the coaxial jet case, the governing equations are made dimensionless using the diameter
$\widetilde{D}$
of the nozzle and the maximum velocity
$\widetilde{U}$
at the inlet boundary of the computational domain
$\unicode[STIX]{x1D6E4}_{in}$
, which is illustrated in figure 2. In order to reproduce the results reported by these authors for
$Re=\widetilde{U}\widetilde{D}/\unicode[STIX]{x1D708}=1000$
, the same axial velocity profile is assumed on
$\unicode[STIX]{x1D6E4}_{in}$
with
$$\begin{eqnarray}\boldsymbol{U}|_{\unicode[STIX]{x1D6E4}_{in}}=U_{z}(r)\hat{\boldsymbol{z}},\quad U_{z}(r)=\tanh \left[5\left(1-r\right)\right],\end{eqnarray}$$
while on the remaining portions of
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{c}$
, the imposed boundary conditions are unchanged with respect to the coaxial jet configuration. The computational domain has the same spatial extent of that adopted by Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013), and contains
$317\,406$
elements, corresponding to
$2\,072\,740$
degrees of freedom. The global mode spectrum obtained for
$m=0$
is illustrated in figure 3(a), and it features all the branches reported in the reference paper. In particular, the branch of modes
$b_{2}$
, which describes the advection of the vortical structures generated by the shear layer instability downstream of the nozzle, is very well reproduced. The computed gains for the optimal perturbations at
$m=0$
are illustrated in figure 3(b) as a function of the time horizon
$\unicode[STIX]{x1D70F}$
. The results are once again in good agreement with the reference ones (reported in the same figure), indicating a monotonic growth of the amplification factor of up to three orders of magnitude.
Figure 3. Modal and non-modal stability analyses of the single jet flow for
$Re=1000$
and
$m=0$
. (a) Global eigenspectrum. (b) Gains associated with the optimal perturbations as a function of the time horizon
$\unicode[STIX]{x1D70F}$
: present results (round circle) compared with those reported by Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013) for the same configuration (continuous line).
4 Base flow
Figure 4. Variation of the base flow velocity profiles through the nozzle for
$R_{u}=1$
and different values of
$Re$
. (a) Axial velocity profile imposed at the inflow and resulting velocity profiles close to the nozzle exit (
$\bar{z}=0.1$
). (b) Radial velocity profiles close to the nozzle exit (
$\bar{z}=0.1$
).
Before addressing the analysis of the linearised perturbation dynamics, the main properties of the steady, axisymmetric base flow produced by the two coaxial jets are shortly described. As mentioned, the base flow represents a numerical approximation to an exact steady solution of the Navier–Stokes equations for the present flow configuration. The inflow profiles are chosen to approximate those measured in experiments employing converging nozzles (Talamelli & Gavarini Reference Talamelli and Gavarini2006) and have a top-hat shape similar to that employed in Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013) for laminar flow analyses. It is worth mentioning that part of the stability analysis has also been performed using fully developed Poiseuille velocity profiles at the inflow and the results, not presented in this work for conciseness, were qualitatively unchanged.
By inspecting the base flow field, it can be observed that the axial velocity profile prescribed on the inlet boundary
$\unicode[STIX]{x1D6E4}_{in,i}\cup \unicode[STIX]{x1D6E4}_{in,o}$
gradually modifies through the nozzle due to the boundary layer growth inside the pipes. A comparison between the axial velocity profiles at the inlet and those resulting at the outlet of the nozzle is illustrated in figure 4 for
$R_{u}=1$
and different values of
$Re$
. At the pipe outlet, the outer jet displays a nearly parabolic profile with a well-defined peak, while for the inner fluid stream, the ‘top-hat’ profile shape is almost preserved. As expected, these shape modifications become more pronounced as the Reynolds number is reduced. For both fluid streams the maximum inlet velocity is increased up to the
$30\,\%$
at the nozzle exit, while the ratio between the maxima is not significantly altered, due to a comparable growth of the boundary layers. The exact values of the outlet–inlet maximum velocity ratios for the inner and outer fluid streams,
$R_{u,e}^{i}$
and
$R_{u,e}^{o}$
respectively, are reported in table 2. Figure 4(b) also reports the radial component of the velocity field at the pipe outlets as a function of the Reynolds number. It can be observed that
$U_{r}$
is about one order of magnitude smaller than
$U_{z}$
and that the modulus of
$U_{r}$
is reduced by increasing the Reynolds number.
Figure 5. Base flow axial velocity field in the proximity of the nozzle for
$Re=1000$
and different values of
$R_{u}$
. White lines are used to represent the flow streamlines while the grey shaded area corresponds to the duct wall separating the two coaxial fluid streams. (a)
$R_{u}=0.5$
. (b)
$R_{u}=1$
. (c)
$R_{u}=1.5$
. (d)
$R_{u}=2$
. The same colour scale is employed for all the panels and the colour bar is reported in figure (d).
Table 2. Outlet–inlet ratios of the maximum axial velocity component characterising the base flow for
$R_{u}=1$
and different values of
$Re$
. See also figure 4 for
$\bar{z}=0.1$
.
The presence of a thick duct wall separating the two coaxial jet streams is responsible for the flow separation at the nozzle exit, with the formation of a closed region of recirculating fluid, as in the case of steady bluff-body wakes. The local structure of the base flow is illustrated in figure 5 for
$Re=1000$
and four different values of
$R_{u}$
. For
$R_{u}=0.5$
the wake is nearly symmetric with respect to the duct wall centreline. However, by increasing
$R_{u}$
the flow topology undergoes a substantial modification. Both vortical rings attached to the wall gradually shrink and the one located on the outer stream side eventually disappears. In addition, a more elongated region of low-speed fluid in the axial direction is observed and a third slender vortex ring is formed away from the wall. It is worthwhile to note that this change in the base flow topology occurs without any bifurcation of the initial state, similarly to what observed by Carini, Giannetti & Auteri (Reference Carini, Giannetti and Auteri2014a
) in the flow past two side-by-side circular cylinders at low Reynolds numbers.
Figure 6. Axial and radial velocity profiles of the computed base flow for
$Re=1000$
and
$R_{u}=1$
at different stations along the
$z$
-axis. (a)
$U_{z}(r,z)$
. (b)
$U_{r}(r,z)$
.
Although the weakly non-parallel assumption can provide a good approximation in the stability analysis of the considered flow, the separation which occurs behind the duct wall introduces important non-parallel effects, especially in the region close to the nozzle. These effects are better illustrated in figure 6, by comparing the axial and radial velocity profiles extracted from the base flow at
$Re=1000$
and
$R_{u}=1$
, at different stations along the
$z$
-axis. As expected, the radial velocity component is very small if compared with the axial one, being three orders of magnitude lower, figure 6(b). However, the influence of the duct wall wake is still visible almost up to
$z\approx 40$
, figure 6(a), and in the region downstream of the nozzle the variation of the velocity profiles
$U_{z}$
and
$U_{r}$
profiles is not as slow as it is usually assumed to be under the weakly non-parallel assumption. These observations further motivate the adoption of a global approach in the stability analysis of the flow.
5 Stability analysis
5.1 Eigenspectrum
Figure 7. Global stability analysis of the coaxial jet flow. Panel (a) depicts the neutral curve associated with the leading axisymmetric mode in the
$Re$
–
$R_{u}$
plane. The minimum critical Reynolds number,
$Re=1356$
for
$R_{u}\approx 1.27$
, is highlighted with a (red) dash-dotted line. The round markers,
$P_{i}$
, highlight the values at which the spectra in figure 8 have been computed. The dashed line represents the neutral curve as a function of the bulk Reynolds number
$Re_{b}$
, reported on the top
$x$
-axis. Panel (b) reports the Strouhal number of the marginally stable mode along the neutral curve, plotted as a function of the critical velocity ratio
$R_{u}$
. The continuous line corresponds to the frequency scaled with the inner jet velocity,
$St_{s}=s\unicode[STIX]{x1D714}_{c}/(2\unicode[STIX]{x03C0}U_{i})$
, while the dashed line shows the same frequency scaled by the bulk velocity,
$St_{s,b}=s\unicode[STIX]{x1D714}_{c}/(2\unicode[STIX]{x03C0}U_{b})$
, reported on the top
$x$
-axis.
Figure 8. Global eigenspectrum computed at four different points along the neutral curve for
$m=0,1,2$
. (a)
$P_{1}$
,
$Re=5000$
and
$R_{u}=0.505$
. (b)
$P_{2}$
,
$Re=1420$
and
$R_{u}=1.0$
. (c)
$P_{3}$
,
$Re=1386$
and
$R_{u}=1.5$
. (d)
$P_{4}$
,
$Re=1559$
and
$R_{u}=2.0$
.
The computation of the global eigenspectrum of the linearised coaxial jet flow for
$m=0$
reveals the presence of a pair of complex-conjugate eigenvalues which become unstable over a wide range of the governing parameters,
$Re$
and
$R_{u}$
. This region of instability is illustrated in figure 7(a) as a grey shaded area delimited by the neutral curve of the considered mode (black line). The neutral curve is characterised by a minimum value of the critical Reynolds number of
$Re=1356$
achieved for
$R_{u}\approx 1.27$
. On its lower branch, for
$R_{u}<1$
, the dependence of
$R_{u}$
at the instability threshold becomes gradually weaker as
$Re$
is increased beyond
$Re=3000$
, with a critical velocity ratio
$R_{u}\approx 0.5$
. Indeed, at low velocity ratios the flow is dominated by the inner jet instability, approaching the behaviour of a single jet, which displays a stable eigenspectrum for all Reynolds numbers (see, e.g., Garnaud et al.
Reference Garnaud, Lesshafft, Schmid and Huerre2013). It appears that for
$R_{u}\lesssim 0.5$
the outer jet is essentially reduced to a co-flow for the inner jet and its influence is not strong enough to alter the dynamics of the inner jet. This does not imply that
$R_{u}\approx 0.5$
is a strict lower minimum for the velocity ratio: the neutral curve may have an horizontal asymptote for any value below
$0.5$
.
It is observed that the critical base flow presents a pair of counter rotating vortex rings, located in the wake of the separating wall (see figures 5 and 11), all along the lower branch of the neutral curve. While these vortex rings are approximately of the same size for
$R_{u}\approx 0.505$
, the outer vortex ring shrinks when moving along the lower branch of the neutral curve and increasing the velocity ratio, until disappearing completely for
$R_{u}\approx 0.90$
and
$Re\approx 1495$
. On the other hand, the topology of the base flow is only slightly modified along the upper branch of the neutral curve, and presents only one recirculating vortex, similarly to what was described for a subcritical
$Re$
(see figure 5). It therefore seems that the onset of instability of the wake of the separating wall is not related to a topological feature of the flow. The non-dimensional critical mode frequency scaled using the wall duct thickness,
$St_{s}=s\unicode[STIX]{x1D714}_{c}/(2\unicode[STIX]{x03C0}U_{i})$
, is illustrated in figure 7(b) as a function of
$R_{u}$
. Except for its behaviour in the neighbourhood of
$R_{u}=0.5$
, the mode frequency increases almost linearly with
$R_{u}$
in analogy with what has been observed by Talamelli & Gavarini (Reference Talamelli and Gavarini2006). To further investigate the similarity to the results by Talamelli & Gavarini (Reference Talamelli and Gavarini2006), figure 7(b) also reports the Strouhal number referred to the bulk velocity of the two jets,
$St_{s,b}=s\unicode[STIX]{x1D714}_{c}/(2\unicode[STIX]{x03C0}U_{b})$
. It can be observed that, when scaled by the mean velocity, the mode frequency is approximately constant for
$R_{u}>1$
, in agreement with figure 11 in Talamelli & Gavarini (Reference Talamelli and Gavarini2006), but
$St_{s,b}$
does not assume a constant value for low velocity ratios, where non-parallelism effects appear to be more relevant.
Figure 8 illustrates the eigenspectra computed at the stations
$P_{i}$
along the neutral curve of figure 7(a). In the same figure, the eigenvalues obtained for
$m=1$
and
$m=2$
are also reported. For these wavenumbers, i.e.
$m=1,2$
, all the eigenvalues have a negative growth rate, thus confirming that the primary instability is axisymmetric. In addition to the leading eigenvalue (black dot lying on the imaginary axis) and in analogy with the single jet spectrum, figure 3(a), different mode branches are observed in the spectrum for
$m=0$
. More precisely, with reference to figure 8(b), the branch
$b_{1}$
is composed by nearly steady modes associated with vortical structures in the fluid region surrounding the jets, which are typical of parallel flows. In contrast, the modes belonging to the branch
$b_{2}$
are localised within the jet shear layers. Finally the branch
$b_{3}$
contains poorly conditioned eigenvalues, which are not exactly recovered when changing the adopted complex shift in the eigenvalue computations. Such branch structure is common to all the panels of figure 8.
The present results are in good agreement with those by Talamelli & Gavarini (Reference Talamelli and Gavarini2006) for low velocity ratios. It is interesting to notice that the analysis by these authors is mainly focused on null or small back flow in the separating wall region. This seems to suggest that the structure of the counter rotating vortex rings has a minor role for low values of
$R_{u}$
, where the instability may originate from a synchronisation between the characteristic frequencies of the two shear layers separated by the duct wall, and therefore be mostly dependent on the velocity ratio alone. On the other hand, Talamelli & Gavarini (Reference Talamelli and Gavarini2006) observe an upper
$R_{u}$
limit for the instability (sensitive to the boundary layer thickness) which is not found in the present global analysis. This discrepancy seems to suggest that the recirculation region, and the non-parallelism of the flow, plays a more important role on the instability for high velocity ratios. Even Talamelli & Gavarini (Reference Talamelli and Gavarini2006) recognise that ‘for velocity ratio
$R_{u}<0.5$
and
$R_{u}>2.6$
an absolute instability occurs only if a back flow is generated in the wake region. For instance this effect may be obtained by aspirating the flow in proximity of the separating wall’. The present analysis shows that aspiration is not necessary, since the recirculation region behind the separating wall naturally provides a back flow, but indicates the necessity for a fully non-parallel base flow to correctly capture the nature of the instability.
It then appears that the global instability that sets in for
$R_{u}\gtrsim 0.5$
is related to the presence of the blunt wall separating the two jets. This is confirmed by a global, linear stability analysis, not reported in the present work for conciseness, showing that the eigenvalue that becomes unstable is largely stable when the blunt wall is substituted by a streamlined one. However, the physical origin of this instability is not entirely clear, as it is not clear for cylinder flows. No discontinuous modifications in the structure of the base flow can, in fact, be observed when crossing the neutral curve. Besides, the same critical global mode is observed even for high velocity ratios, where the topology of the recirculating region is fundamentally altered, presenting one vortex ring instead of the two observed for low
$R_{u}$
. Moreover, the adopted scaling, where the Reynolds number is insensitive to the velocity of the outer jet, being solely based on the velocity of the inner jet, leads to a slow change of the critical
$Re$
with
$R_{u}$
in the upper branch of the neutral curve. This is not the case if a different scaling is used as, for example, by employing a Reynolds number based on the total bulk velocity of the two jets. In this case, shown by the dashed line in figure 7(a), the dependence of the upper branch of the neutral curve on the velocity ratio is stronger.
5.2 Leading global modes and the wavemaker
Figure 9. Direct global mode for
$m=0$
. Panels (a) and (b) represent the axial and radial velocity components, respectively, of the leading global mode at criticality for
$Re=1420,R_{u}=1.0$
(point
$P_{2}$
in figure 7). The same quantities are depicted in panels (c) and (d) for the critical mode at
$Re=1559,R_{u}=2.0$
(point
$P_{4}$
in figure 7). Note that the figures only reproduce a portion of the computational domain which is always characterised by
$z_{max}\geqslant 40$
. Mesh
$M_{60}$
, with
$z_{max}=60$
, was used for these particular fields. The complete fields can be found at https://doi.org/10.1017/jfm.2017.290.
Figure 10. Critical adjoint global mode for
$m=0$
represented by the magnitude of the associated velocity field. (a)
$Re=1420$
,
$R_{u}=1.0$
. (b)
$Re=1559$
,
$R_{u}=2.0$
. Note that the figures only reproduce a portion of the computational domain which is always characterised by
$z_{min}\leqslant -5$
. Mesh
$M_{60}$
, with
$z_{min}=-5$
, was used for these particular fields. The complete fields can be found at https://doi.org/10.1017/jfm.2017.290.
Figure 11. Structural sensitivity map
${\Vert \unicode[STIX]{x1D61A}_{0}(r,z)\Vert }_{F}$
associated with the critical global mode at four different points along the neutral curve. White lines are used to represent the base flow streamlines. (a)
$P_{1}$
,
$Re=5000$
and
$R_{u}=0.505$
. (b)
$P_{2}$
,
$Re=1420$
and
$R_{u}=1.0$
. (c)
$P_{3}$
,
$Re=1386$
and
$R_{u}=1.5$
. (d)
$P_{4}$
,
$Re=1559$
and
$R_{u}=2.0$
. (e)
$P_{4}$
,
$Re=1559$
and
$R_{u}=2.0$
(extended view). Note that the figures only reproduce a portion of the computational domain which is always characterised by
$z_{max}\geqslant 40$
and
$z_{min}\leqslant -5$
. Mesh
$M_{60}$
, with
$z_{max}=60$
and
$z_{min}=-5$
, was used for these particular fields. The complete fields can be found at https://doi.org/10.1017/jfm.2017.290.
The direct global mode computed for
$R_{u}=1$
is represented in figures 9(a) and 9(b) by means of the real part of its axial and radial velocity components, respectively. This mode is formed by an array of counter-rotating vortex rings developing in the wake of the separating duct wall. The amplitude of these structures grows downstream of the nozzle, in the axial direction, reaching a maximum value at
$z\approx 11$
, after which they slowly decay. The spatial structure of this eigenmode resembles that of the critical eigenmode in the wake of a circular cylinder (Sipp & Lebedev Reference Sipp and Lebedev2007). Direct numerical simulations, in fact, show an ‘annular’ alternating vortex street originating in the wake of the duct wall. Such spatial structure substantially differs from that characterising the shear layer modes belonging to branch
$b_{2}$
, which display an exponential spatial growth up to the outlet boundary of the computational domain, similarly to what observed by Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013) for a single jet. The structure of the critical mode substantially changes when moving along its neutral curve. Figure 9(c,d) represent the critical mode for
$R_{u}=2,Re=1559$
. The maximum of the amplitude is located further downstream, at
$z\approx 20$
, and the shape is qualitatively modified. This modification is reflected in the DNS which shows a Kelvin–Helmholtz instability in the shear layer between the jets and at the interface between the outer jet and the quiescent surrounding fluid.
Differently from the direct mode, the adjoint mode, illustrated in figure 10, results essentially concentrated within the nozzle. Its spatial distribution is found in the form of upstream inclined structures localised inside the pipe boundary layers, achieving a maximum intensity close to the sharp corners of the duct wall. When the velocity ratio is changed, the adjoint mode moves from the inner pipe, where it is mainly located for low
$R_{u}$
, to the outer pipe, for
$R_{u}>1$
.
The knowledge of both the direct and the adjoint modes allows us to identify the regions of the flow where the considered instability originates, according to the wavemaker analysis introduced by Giannetti & Luchini (Reference Giannetti and Luchini2007). The ‘core’ of the global instability can indeed be associated with the maximum values of the function
$\Vert \unicode[STIX]{x1D61A}_{0}(r,z)\Vert _{F}$
, defined in § 2. Figure 11 illustrates the sensitivity map computed along the neutral curve. Except for the case at
$R_{u}=2$
, the spatial distribution of
$\Vert \unicode[STIX]{x1D61A}_{0}(r,z)\Vert _{F}$
results highly localised behind the duct wall, featuring negligible values both inside the pipes and far from the nozzle. For
$R_{u}=1$
, figure 11(b), the maximum intensity is attained within a ‘double lobe’ structure which is nearly symmetric with respect to the duct wall centreline. A similar shape of the wavemaker has been observed for several two-dimensional bluff-body flow configurations at low Reynolds numbers (Giannetti & Luchini Reference Giannetti and Luchini2007; Pralits et al.
Reference Pralits, Brandt and Giannetti2010; Carini et al.
Reference Carini, Giannetti and Auteri2014b
). Either decreasing or increasing the velocity ratio, the wavemaker structure becomes distinctly asymmetric, as shown in figure 11(a,c,d), featuring a single pocket of maximum sensitivity located on the side of the faster flowing jet. Finally, by increasing the velocity ratio up to
$R_{u}=2$
, another peculiar modification of the sensitivity map occurs. As already mentioned, in this case small but non-negligible values of the norm of the sensitivity tensor are observed over a wide region downstream of the nozzle. A similar behaviour has also been reported by Carini et al. (Reference Carini, Giannetti and Auteri2014a
) in a preliminary analysis of the secondary bifurcation of the flow past two side-by-side circular cylinders. However, differently from that case, in the present configuration the region of non-zero sensitivity values is entirely contained within the adopted computational domain. An extended view is presented in figure 11(e), and ensures the proper convergence of the computed mode with respect to the downstream location of the outlet boundary. This is also confirmed by the convergence study reported in the appendix A. Although the wavemaker extends downstream of the nozzles for high velocity ratios, the most significant contribution to the instability remains localised close to the separating wall. This observation is confirmed by the work of Tammisola (Reference Tammisola2012) via a global and local analysis. The author reports that, for a pair of two-dimensional jets, the downstream portion of the wavemaker is not detected by a local stability analysis. This indicates that the ‘true’ wavemaker is the one close to the pipe exit, while the downstream area of sensitivity is associated with the upstream reflection of perturbations generated by the ‘true’ wavemaker. In the present analysis the magnitude of the downstream portion of the wavemaker is lower than that of the upstream portion, indicating that this effect is not as prominent as in Tammisola (Reference Tammisola2012). The resemblance between the wavemakers in the two works, however, suggests that even for coaxial jets at high velocity ratios the downstream portion of the wake does not have a merely convective function, but is coupled to the near wake dynamics.
5.3 Transient growth
Figure 12. Gains associated with the optimal perturbations as a function of
$\unicode[STIX]{x1D70F}$
for
$m=0$
and different values of
$R_{u}$
. (a)
$Re=1000$
. (b)
$Re=1250$
. The dashed (blue) line shows the amplification of the same optimal perturbations attained on the nonlinear flow
$R_{u}=1$
.
Figure 13. Spatial structures associated with the evolution of the optimal perturbations at
$Re=1000$
depicted by means of the kinetic energy distribution
$|\boldsymbol{u}|^{2}/2$
for different velocity ratios: (a,b)
$R_{u}=0.5$
, (c,d)
$R_{u}=1.0$
and (e,f)
$R_{u}=1.5$
. In (a,c,e) the optimal initial condition is displayed while in the (b,d,f) the resulting linearised flow response at the maximum amplification time (black filled symbols in figure 13
a) is illustrated. Note that the figures only reproduce a portion of the computational domain which is always characterised by
$z_{max}\geqslant 40$
and
$z_{min}\leqslant -5$
. Meshes
$M_{60}$
and
$M_{80}$
, with
$z_{min}=-5$
and
$z_{max}=60$
and
$80$
respectively, were used for these particular fields.
The linearised flow dynamics is now investigated by means of a non-modal stability analysis. The non-modal analysis is here restricted to axisymmetric perturbations in the subcritical regime. Figure 12 illustrates the optimal gain curve
$G_{0}^{opt}(\unicode[STIX]{x1D70F})$
computed for three different values of the velocity ratio, at
$Re=1000$
and
$Re=1250$
(figures 12(a) and 12(b), respectively). For
$R_{u}=0.5$
the amplification factors are of the same order as those characterising the single jet, already shown in figure 3(b). For low velocity ratios the transient amplification of the linearised perturbations is essentially driven by the inner shear layer response. In contrast, huge gains are obtained for a unitary velocity ratio and
$G_{0}^{opt}(\unicode[STIX]{x1D70F})$
further increases up to
$20$
orders of magnitude for
$R_{u}=1.5$
. When
$R_{u}=1.0$
, the maximum amplification is reached at
$\unicode[STIX]{x1D70F}=60$
, which approximately corresponds to
$z=34$
, while for
$R_{u}=1.5$
, optimal perturbations are still being amplified when they are advected outside of the computational domain.
These results are better illustrated in figure 13 for each considered velocity ratio and
$Re=1000$
. This figure depicts both the optimal initial perturbation and the resulting linearised flow field at the maximum amplification time (corresponding to the black filled symbols in figure 12
a). In analogy with the critical adjoint mode shape, the optimal perturbations are found to be highly localised at the corners of the duct wall. In contrast to the adjoint modes, though, the regions of maximum intensity are observed to be localised on the side of the stream with lower maximum velocity for
$R_{u}\neq 1$
, while they are located on both jet sides for a unitary velocity ratio. The difference between the results of the wavemaker analysis and those of the transient-growth analysis might seem ‘surprising’ at first: the wavemaker displays higher sensitivity to perturbations on the side of the faster flowing jet, while the optimal initial perturbations are located on the inner side of the duct walls. It should be noted, however, that the two analyses correspond to two different kinds of instability. The structural sensitivity is associated with an intrinsic instability mechanism: the base flow is unstable independently of external perturbations; the transient growth, instead, is related to the input/output relationship between the perturbations introduced into the flow and how they are transformed by the flow itself. This fundamental difference is reflected in the results of the present analysis: the absolute instability is associated with the presence of a blunt separating wall and to the interaction of the shear layers separating from it. In contrast, the convective instability is associated with the transient growth of perturbations seeded within the inlet ducts, near their walls, and amplified in the entrance region by a Orr mechanism (Garnaud et al.
Reference Garnaud, Lesshafft, Schmid and Huerre2013; Sipp & Marquet Reference Sipp and Marquet2013) and in the jet region by the Kelvin–Helmholtz instability mechanism. It is interesting to notice, here, the more effective amplification of perturbations seeded near the external wall of the jets.
To assess the relevance of the linear optimal perturbations on the actual flow, the same perturbations have been employed as initial condition in nonlinear simulations. It should be noted that this operation is not connected to a nonlinear optimals analysis; it is simply an exercise to assess the actual amplification of these perturbations. It should not be forgotten that the nonlinear flow depends on the initial amplitude. The results are illustrated in figure 12(a) for
$R_{u}=1.0$
and
$Re=1000$
setting the initial amplitude of the velocity perturbation field to a value between 0.4 % and 2 % of the base flow velocity, measured in the
$L^{\infty }$
norm of the axial component. As expected, the gains are lower in the nonlinear flow, and the kinetic energy saturates for
$\unicode[STIX]{x1D70F}>30$
, forming a plateau. The amplification of these initial conditions, though, is still considerably high, reaching a maximum value of
$1.3\times 10^{5}$
, two orders of magnitude higher than the linear maximum for a single round jet. This is especially true since quite large initial perturbations have been considered. Despite these enormous growths, the nonlinear flow confirms its absolute stability, and returns to the steady state after advecting the perturbations outside of the domain.
Figure 14. Flow snapshot of the fully developed unsteady coaxial jet flow obtained from three-dimensional DNS at
$R_{u}=1$
and
$Re=1500$
: iso-surfaces of the
$\unicode[STIX]{x1D706}_{2}$
criterion (
$\unicode[STIX]{x1D706}_{2}=-0.5$
) coloured by the azimuthal component of the vorticity field. A planar cut of the azimuthal vorticity field is also illustrated. The light grey shaded structure corresponds to the nozzle geometry.
6 Direct numerical simulations
Figure 15. Time history and frequency content of the radial velocity component extracted from the DNS of the axisymmetric coaxial jet flow for
$R_{u}=1$
and two different Reynolds numbers, i.e.
$Re=1000$
(a,d,g) and
$Re=1500$
(b,e,h). The velocity probes are located at
$z=2.0$
and
$r=0.25,0.55,0.8$
corresponding to (a–c), (d–f) and (g–i), respectively. The third column shows the spectral content of the velocity signals for the case at
$Re=1500$
, for which a periodic solution is established.
Direct numerical simulations of the two coaxial jets have been performed in order to assess the onset of the unstable axisymmetric mode predicted by the linear global stability analysis. Since the corresponding instability mechanism is expected to be dominant only in the range of unitary velocity ratios, the DNS computations have been focused on the value of
$R_{u}=1$
. Figure 14 illustrates the vortex structures developing in the unsteady fully three-dimensional flow field at
$Re=1500$
. The structures are represented by iso-surfaces of the
$\unicode[STIX]{x1D706}_{2}$
criterion, for
$\unicode[STIX]{x1D706}_{2}=-0.5$
, and coloured with the corresponding value of the azimuthal component of the vorticity field. As can be observed in this figure, the flow does not feature any helical mode, thus confirming the expected axisymmetric nature of the leading instability. Based on these results, further DNS investigations are carried out under the axisymmetric constraint, which allows us to considerably reduce the computational cost. The simulations have been impulsively started from quiescent conditions and then advanced in time with a non-dimensional time step of
$\unicode[STIX]{x0394}t=1\times 10^{-3}$
up to the achievement of either a steady or unsteady condition. In order to evaluate the frequency of the instability, virtual probes measuring the velocity components are introduced in the simulation. The measurements are presented in figure 15 for three of the probes located at
$z=2.0$
and different radial positions:
$r=0.25,0.55$
and
$0.8$
, corresponding approximately to the inner jet, the duct wall wake and the external jet. For
$Re=1000$
(a,d,g), all the unsteady perturbations exponentially decay and the flow converges towards the steady base state, in agreement with the instability threshold of
$Re\approx 1420$
defined by the neutral curve of figure 7(a). In contrast beyond the neutral curve, for
$Re=1500$
(figure 15
b,e,h) a stable periodic solution is established. The frequency content associated with the probe signals (figure 15
c,f,i) is characterised by a main peak occurring at the same non-dimensional frequency of
$St_{s}=0.0924$
for all the probes. Such value of the vortex shedding frequency is in very good agreement with the unstable global mode frequency of
$St_{s}=0.0912$
resulting from the linear stability analysis, thus supporting the proposed interpretation of the ‘global’ nature and origin of the considered flow instability. The higher frequency peaks appearing in the spectra of figure 15(c,f,i) correspond to second- and third-order harmonics with
$2St_{s}=0.1847$
and
$3St_{s}=0.2783$
.
7 Summary and conclusions
The linear stability of the incompressible flow produced by two coaxial fluid streams separated by a thick duct wall has been investigated within a fully non-parallel framework in order to properly account for the effects of the solid wall and of the jet spreading. For such a configuration, and differently from the case of a single round jet, the flow becomes globally unstable due to an oscillatory axisymmetric mode associated with the vortex shedding in the wake of the duct wall. This finding agrees with the existence of a finite region of absolute instability predicted by the local stability analysis (Talamelli & Gavarini Reference Talamelli and Gavarini2006). Despite the low Reynolds numbers that characterise it, the instability mechanism discovered in the present analysis is known to dominate the flow for nearly unitary velocity ratios for higher Reynolds numbers (see, e.g. Segalini & Talamelli Reference Segalini and Talamelli2011). The characteristic annular alternating vortex street is, in fact, dominating the large-scale motion up to at least
$Re=13\,800$
, while being modulated by small-scale turbulence. It was found that the region of linear instability corresponding to the global mode extends over a wide range of the governing parameters, with a minimum critical Reynolds number of
$Re\approx 1356$
(
$R_{u}\approx 1.27$
) and a lower threshold of
$R_{u}\approx 0.5$
, which is almost independent of
$Re$
. The spatial structure of the leading mode results remarkably different from that characterising the jet shear layer modes, the former developing upstream of the latter and featuring a well-defined maximum within the computational domain. By inspecting the direct-adjoint product, the ‘core’ of this instability results located in the separated flow region which forms in the near wake of the duct wall, showing two pockets of maximum sensitivity on the two sides of the wake of the wall. As the velocity ratio is increased or decreased, the shape of the wavemaker is modified, featuring only one ‘lobe’ located on the side of the faster flowing stream.
A three-dimensional DNS at
$R_{u}=1$
and
$Re=1500$
has been performed to assess the axisymmetric nature of the unsteady flow and further axisymmetric simulations have been used to compare the vortex shedding frequency with the global mode frequency. Good agreement was found for all values of the parameters, with a slight difference of approximately
$1.3\,\%$
.
Finally an optimal perturbation analysis has been carried out to complete the linear global stability description of the flow. Quite surprisingly, huge gain factors are found to characterise the transient response of the linearised flow dynamics at high velocity ratios. The optimal growth far exceeds that of a single round jet, achieving an amplification of up to
$20$
orders of magnitude. In analogy with the results of Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013), optimal perturbations result highly localised around the sharp corners of the nozzle. In particular it is interesting to note that, as the velocity ratio is varied and differently from the wavemaker, their location shifts towards the side of the lower-speed fluid stream. The same linear optimal perturbations were also used as initial conditions in nonlinear simulations. In spite of the expected reduction in the amplification attained by these perturbations, the gains remain considerably high, reaching values which are two orders of magnitude higher than the linear maximum for a single round jet.
Figure 16. Branch
$b_{2}$
of the global eigenspectrum computed using meshes characterised by a different streamwise extent of the domain.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2017.290.
Appendix A. Sensitivity to domain size and grid resolution
In order to investigate the influence of the domain size on the stability results, several computations have been performed by varying both the radial and the axial extension of the domain. All computations were carried out with the same local mesh size,
$\unicode[STIX]{x0394}\ell$
, varying from
$\unicode[STIX]{x0394}\ell =0.4$
close to
$\unicode[STIX]{x1D6E4}_{o}$
, down to a minimum value of
$\unicode[STIX]{x0394}\ell =0.009$
within the refinement zone (blue shaded area in figure 1), these sizes of the elements have been determined following a mesh refinement independence study; the different features of the employed meshes are listed in table 1. As occurs in the case of the single jet (Garnaud et al.
Reference Garnaud, Lesshafft, Schmid and Huerre2013), the stable eigenvalues belonging to the branch
$b_{2}$
of the spectrum, which are associated with exponential downstream growing modes, are strongly influenced by the position of the outlet boundary. This is confirmed in figure 16, which reports the branch
$b_{2}$
when computed by using meshes characterised by an increasing axial length. When
$z_{max}\geqslant 100$
(meshes
$M_{100}$
–
$M_{160}$
) the overall displacement of this branch becomes smaller but the eigenvalues that constitute it are joined by other ill-conditioned eigenvalues. In contrast, the leading unstable eigenvalue displays a good convergence with respect to the spatial extension of the domain: the first four significant digits of its frequency are, in fact, independent on the mesh length, as shown in table 3 for
$R_{u}=1$
and
$Re=1420$
. Similarly, table 4 lists the variation of the amplification factor
$G_{0}^{opt}(\unicode[STIX]{x1D70F})$
computed using the different grids for
$Re=1000$
and three different values of
$\unicode[STIX]{x1D70F}$
. It is possible to observe that
$G_{0}^{opt}(\unicode[STIX]{x1D70F})$
shows very little sensitivity to domain size. This is explained by the fact that all optimal perturbations are located in the nozzle region, with no contributions from other areas of the domain.
Table 3. Computed leading eigenvalue for
$R_{u}=1$
and
$Re=1420$
using computational domains of different length.
Table 4. Variation of the computed gain factor
$G_{0}^{opt}(\unicode[STIX]{x1D70F})$
for
$R_{u}=1$
and
$Re=1000$
obtained using different computational domains, with respect to the reference
$M_{100,z_{i}15}$
grid. The variation is also normalised with respect to the values obtained with
$M_{100,z_{i}15}$
.
