Artificial eigenmodes in truncated flow domains

TL;DR

This paper investigates how artificial boundary conditions in truncated flow domains generate spurious eigenmodes, affecting the interpretation of open flow stability analyses, and proposes methods to mitigate these artifacts.

Contribution

It identifies the origin of spurious eigenmodes caused by boundary feedback in truncated domains and demonstrates how absorbing layers can reduce these artifacts.

Findings

Spurious eigenmodes depend on domain length and boundary conditions.

Pressure feedback from outflow to inflow decays with domain length cubed.

Absorbing layers near outflow reduce spurious feedback.

Abstract

Whenever linear eigenmodes of open flows are computed on a numerical domain that is truncated in the streamwise direction, artificial boundary conditions may give rise to spurious pressure signals that are capable of providing unwanted perturbation feedback to upstream locations. The manifestation of such feedback in the eigenmode spectrum is analysed here for two simple configurations. First, explicitly prescribed feedback in a Ginzburg-Landau model is shown to produce a spurious eigenmode branch, named the 'arc branch', that strongly resembles a characteristic family of eigenmodes typically present in open shear flow calculations. Second, corresponding mode branches in the global spectrum of an incompressible parallel jet in a truncated domain are examined. It is demonstrated that these eigenmodes of the numerical model depend on the presence of spurious forcing of a local k+ instability wave at the inflow, caused by pressure signals that appear to be generated at the outflow. Multiple local k+ branches result in multiple global eigenmode branches. For the particular boundary treatment chosen here, the strength of the pressure feedback from the outflow towards the inflow boundary is found to decay with the cube of the numerical domain length. It is concluded that arc-branch eigenmodes are artifacts of domain truncation, with limited value for physical analysis. It is demonstrated, for the example of a non-parallel jet, how spurious feedback may be reduced by an absorbing layer near the outflow boundary.

Figures

Tables

Table 1. Table 1: Numerical configurations of the various eigenmode calculations discussed in this section.

	$x_{min}$	$x_{max}$	upstream BC	absorbing layer
case 1	0	40	Dirichlet	no A
case 2	$-5$	40	Dirichlet	no
case 3	$-10$	40	Dirichlet	no
case 4	0	40	stress-free	no
case 5	$-5$	40	stress-free	no
case 6	$-10$	40	stress-free	no
case 7	$-5$	40	Dirichlet	yes, ${\lambda }_{max}=0.5$
case 8	$-5$	50	Dirichlet	yes, ${\lambda }_{max}=0.75$

Equations

$${\partial }_t\psi =-U{\partial }_x\psi +\mu (x)\psi +\gamma {\partial }_{xx}\psi +f(x,\psi ).$$

$$\mu (x)=\frac{U^2}{8}\left(1-\frac{x}{20}\right).$$

$$f(x,\psi )=C\exp \left(-\frac{(x-x_a{)}^2}{0.1^2}\right)\psi (x_s),$$

$${\omega }_n=i\left\{\frac{U^2}{8}-\frac{U^2}{4\gamma }+{\gamma }^{\frac{1}{3}}\frac{U^{\frac{4}{3}}}{160^{\frac{2}{3}}}{\zeta }_n\right\},$$

$${\partial }_t\psi =\mathbf{L}\psi \to \Vert (\text{i}\omega \mathbf{I}+\mathbf{L}{)}^{-1}\Vert ={ϵ}^{-1},$$

$$0$$

$${\partial }_t{u}_r$$

$${\partial }_t{u}_x$$

$$u_r=u_x=0\text{at }x=0,$$

$$R{e}^{-1}{\partial }_x{u}_x-p=0\text{at }x=x_{max},$$

$$u_r={\partial }_r{u}_x=0\text{at }r=0\text{ and }r=r_{max}.$$

$$U(r)=\frac{1}{2}\left(1+\tanh \left[\frac{1}{4\theta }\left(\frac{1}{r}-r\right)\right]\right),$$

$$\Delta \hat{p}=-2{\partial }_{r}U{\partial }_x{\hat{u}}_r,$$

$${\hat{p}}_j=\left[A_j{I}_0\left(\frac{j\pi }{L}r\right)+B_j{K}_0\left(\frac{j\pi }{L}r\right)\right]\sin \left(\frac{j\pi }{L}x\right),j\in \mathbb{N},$$

$$-\text{i}\omega \hat{q}=C\hat{q}{\text{e}}^{\text{i}kL}\to C=-\text{i}\omega {\text{e}}^{-\text{i}kL},$$

$$U(r)=\tanh \left[\frac{5}{2}\left(\frac{1}{r}-r\right)\right]$$

$$-i\omega \mathbf{B}\hat{q}=\mathbf{L}\hat{q},$$

$$-i\omega \mathbf{B}\hat{q}=\left[\mathbf{L}-\lambda (x)\mathbf{B}\right]\hat{q}.$$

$$\lambda(x)=\left\{\begin{array}[]{rl}0&\text{~{}~{}~{}for~{}~{}}x\leqslant 19,\\ 0.00625x^{2}-0.2375x+2.25625&\text{~{}~{}~{}for~{}~{}}19<x<21,\\ 0.025(x-20)&\text{~{}~{}~{}for~{}~{}}x\geqslant 21.\end{array}\right.$$

$$A(x)={\left[{\int }_0^1\left(|{\hat{u}}_r{|}^2+|{\hat{u}}_x{|}^2\right)r\text{d}r\right]}^{\frac{1}{2}},$$

Full text

1 Introduction

Linear instability analysis of open flows today is commonly carried out in a so-called ‘global’ framework, where at least two non-homogeneous spatial directions of a steady base state are numerically resolved. In contrast to ‘local’ analysis, where the base state is assumed to be parallel, and unbounded in the flow direction, a global discretisation of an open flow problem in a truncated numerical domain necessitates the formulation of artificial streamwise boundary conditions for flow perturbations. The question then arises in how far such boundary conditions influence the instability dynamics of the truncated flow system.

This paper investigates the effect of spurious pressure feedback, due to domain truncation, on the eigenmode spectrum of incompressible open flow problems. The investigation is motivated by observations made in recent linear instability studies of jet flows Garnaud:2013p1182 ; Coenen2017 , where a prominent family of eigenmodes (black symbols in Fig. 1), referred to as the ‘arc branch’ from here on, was found to present features that suggest a resonance between the inflow and outflow boundaries. Such branches are in fact ubiquitous in many, if not all, global spectra of truncated open shear flows found in the literature: boundary layers (Ehrenstein:2005p1026, ; Ehrenstein:2008p1126, ; Akervik:2008p1176, ), cylinder wakes (Sipp2007, ; Marquet:2008p1127, ), jets (Nichols:2011p1075, ; Garnaud:2013p1182, ), plumes (ChakriPhD, ), three-dimensional boundary layers with roughness elements (Loiseau2014investigation, ; Kurz2016, ) — all these and many others present similar characteristic branches of eigenvalues that are often described as being highly dependent on the type or position of outflow boundary conditions. Typically, no convergence with respect to the length of the numerical box can be attained for such modes.

Ehrenstein & Gallaire Ehrenstein:2005p1026 remark on the resemblance between arc-branch-type global eigenmode structures obtained for a flat-plate boundary layer and spatial modes as found in a local analysis. Åkervik et al. Akervik:2008p1176 document the dependence of boundary layer eigenvalues on the type of outflow boundary conditions. Heaton et al. heaton characterise arc branch modes in the spectrum of a Batchelor vortex as artifacts, speculating that these arise from the limited precision of their numerical scheme, not from domain truncation. Cerqueira & Sipp cerqueira demonstrate that such precision errors indeed lead to spurious ‘quasi-eigenmodes’, but that these differ from the arc branch. In their analysis, quasi-eigenmodes appear in regions of the complex frequency plane where pseudospectrum $\epsilon$ values TrefethenBook are very small, below approximately $10^{-12}$. These modes, also visible in the lower left corner of Fig. 1, are very sensitive to the numerical scheme, to mesh refinement and to the eigenvalue shift parameter. Arc branch modes are found to be robust with respect to those detailscerqueira , but strongly dependent on the numerical domain length. Coenen et al.Coenen2017 show that arc branch eigenfunctions of a jet are characterised by an integer number of wavelengths between the inflow and outflow boundaries, and they suggest an analogy with acoustic modes in a pipe of finite length. This analogy implies that arc branch modes are the result of unwanted resonance between the numerical boundaries, potentially leading to a spurious instability of the numerical system. The present study expands on all these observations, and it aims at a detailed characterisation of unphysical resonance due to imperfect boundary conditions in open shear flow calculations.

The possibility that global instability in truncated systems may be brought about by spurious pressure feedback from boundary conditions was probably first described by Buell & HuerreBuellHuerre . In direct numerical simulations of perturbations in a mixing layer, unstable perturbation growth was observed at long times, in a configuration that was only convectively unstable in a local sense. Such behaviour is inconsistent with the interpretation of local convective instability, and it was demonstrated to be caused by unphysical pressure perturbations emanating from the outflow boundary, which in turn provoked the formation of vortical perturbations at the inflow boundary. Chomaz (Chomaz2005, , section 3.2.2) interprets pressure feedback as a non-local operator variation, arguing that strong non-normality in the spectrum of convection-dominated flows is likely to induce a high sensitivity of global eigenmodes with respect to such feedback.

While the spurious generation of acoustic pressure waves from artificial boundary conditions is an important and much-discussed problem in compressible flow simulations, especially those that aim to accurately capture the acoustic radiation from shear flows (C04, ), the question of how such artifacts may affect the global stability behaviour remains largely unexplored, both in compressible and in incompressible configurations.

The problem is approached here in the following manner. Section 2 presents a global instability analysis of the Ginzburg–Landau equation with explicit feedback from a downstream sensor to an upstream actuator, in order to examine the effect of such feedback on the eigenvalue spectrum in a controlled setting. This model study fully describes the suspected mechanism behind arc branch modes. In Sect. 3 the analysis is extended to a parallel jet flow in a finite-size numerical domain. Global eigenmodes are projected onto their local counterparts, the global pressure field is examined, and spurious feedback effects are analysed. The question how spurious feedback may be reduced in a practical manner is addressed in Sect. 4, for the example of a spatially developing jet.

2 A Ginzburg–Landau model problem

The hypothesis put forward by Coenen et al.Coenen2017 , that the arc branch may be the manifestation of a non-physical upstream scattering of perturbations from the outflow boundary, is first investigated with the help of a simplified model. The linear Ginzburg–Landau equation is written as

[TABLE]

The complex scalar variable $\psi(x,t)$ is a function of time $t$ and of one single spatial coordinate $x$. Constant parameters $U=6$ and $\gamma=1-\text{i}$ are chosen ad hoc, whereas the coefficient $\mu$ varies linearly in $x$ as

[TABLE]

With this particular variation of $\mu$, the system is marginally absolutely unstable at $x=0$, convectively unstable for $0<x<20$, and stable for $x>20$, all in a local sense. Ginzburg–Landau systems of this form, with linearly decreasing $\mu(x)$ and without feedback, $f\equiv 0$, have been extensively used to model the global instability behaviour of spatially developing flows in semi-infinite domains (see for instance Couairon97, ). The instability characteristics mimic, in a very simple and qualitative manner, those of a spreading jet. However, the system (1) is not intended here to predict or reproduce the dynamics of any specific flow.

A forcing term $f(x,\psi)$ is added in (1) in order to provide an explicit closed-loop forcing between the upstream and the downstream end of the flow domain. Some aspects of closed-loop forcing in the Ginzburg–Landau equation are discussed by Chomaz Chomaz2005 ; in the present context, it is used to model a suspected spurious feedback in global shear flow computations. Taken to be of the form

[TABLE]

a forcing proportional to $\psi(x_{s})$ is applied in a close vicinity of $x_{a}$, such that a feedback loop is established between a (downstream) sensor location $x_{s}$ and an (upstream) actuator location $x_{a}$. The complex coefficient $C$ governs the amplitude and phase of the feedback, and the Gaussian spreading around $x_{a}$ is introduced for reasons of numerical resolution.

Equation (1) is discretised on an interval $0\leqslant x\leqslant 40$, with a step size $\Delta x=0.1$, using an upwind-biased seven-point finite difference stencil for the spatial derivatives. A homogeneous Dirichlet boundary condition for $\psi$ is imposed at the upstream boundary, consistent with typical jet conditions. Actuator and sensor locations are chosen close to the boundaries, at $x_{a}=1$ and $x_{s}=39$, where spatial derivatives are well resolved.

Temporal eigenmodes of (1) are sought in the form $\psi(x,t)=\hat{\psi}(x)\,\text{e}^{-i\omega t}$. The eigenvalues of the system without feedback, $C=0$, are known analytically to be

[TABLE]

where $\zeta_{n}$ is the $n^{th}$ root of the Airy function CHR88 . These values are represented in Fig. 2 as circles. Each frame 2(a–d) also displays the numerically computed eigenvalues of systems with feedback, shown as bullet symbols, for different non-zero values of the coefficient $C$. Already with very low-level feedback, $C=10^{-10}$, the spectrum is clearly affected: only the leading three eigenvalues of the unforced system are recovered, and the lower part of the spectrum is replaced with two new stable branches. These branches, named feedback branches in the following, move upward in the complex $\omega$ plane as the feedback coefficient is increased, masking more and more of the original unforced eigenvalues. Note that those affected original eigenvalues are not merely altered by the feedback, but they rather disappear abruptly from the spectrum as they fall below the new branches.

The spectra obtained with $C=10^{-6}$ and $10^{-2}$ resemble those of pure helium jets, as shown in figure 6 of Coenen2017 , and those of cylinder wakes obtained by Marquet et al.(Marquet:2008p1127, , their figure 16). The least stable original eigenvalue lies just above the feedback branches, apparently unaffected. In the strongest feedback case, $C=10^{2}$, the feedback branches have merged into one, overarching the entire spectrum of the feedback-free system, which has altogether disappeared. This picture (Fig. 2$d$) resembles the spectrum of the slowly developing jet of Coenen et al.Coenen2017 , reproduced in our Fig. 1.

The similarity between the feedback branches in the present model and the arc branch in the jet spectrum (Fig. 1) is also manifest in the pseudospectra. If feedback is thought of as a variation of the operator Chomaz2005 , it should be possible to relate feedback-induced eigenvalues to the pseudospectrum. The pseudospectrum is defined here by the spectral norm of the resolvent operator,

[TABLE]

at any complex frequency $\omega$ TrefethenBook . Pseudospectrum $\epsilon$-contours of the Ginzburg–Landau equation are compared between the unforced setting with $C=0$, Fig. 3($a$), and the forced setting with $C=10^{-6}$, Fig. 3($b$). It is seen that the feedback branches align closely with an isocontour of the feedback-free pseudospectrum. This criterion also applies to all other cases displayed in Fig. 2, with different values of $C$, and it is fully consistent with observations made in three flow configurations heaton ; cerqueira ; Coenen2017 . Furthermore, the pseudospectrum of the system with feedback is identical to that of the system without feedback above the feedback branches, whereas the pseudospectrum is flattened, nearly constant, below the feedback branches. The same behaviour is found in the jet pseudospectrum shown in Fig. 1.

Finally, the discrete distribution of feedback modes along the branch is investigated. The strong feedback case $C=10^{2}$ is considered for illustration. Eigenfunctions of successive feedback modes are presented in figure 4, analogous to the representation of jet results in figure 5 of Coenen2017 . The absolute value of $\psi$ is traced in logarithmic scale as a function of $x$. The phase of $\psi$ is always chosen such that the real part of $\psi$ is zero in the sensor location $x_{s}=39$, with the exception of the first mode (Fig. 4$a$), which has no interior wave nodes. Only positive real frequencies are considered.

Clearly, the eigenfunctions are characterised by integer numbers of wavelengths between the actuator and sensor locations, starting from zero wavelengths in the case of the lowest real frequency (Fig. 4$a$, only amplitude variations but approximately constant phase), to one wavelength in Fig. 4$(b)$, and so forth with a continuously increasing count. The series could be continued further along the entire arc branch. The present model seems to reproduce very well the characteristics of arc branch eigenfunctions, as displayed in figure 5 of Coenen2017 .

The position and spacing of eigenvalues along the branch curve, fixed approximately by an isocontour of the pseudospectrum, appears to be determined by a fitting phase relation between $\psi$ at the actuator location and the applied feedback. To further illustrate this mechanism, the phase of the feedback coefficient $C$ is varied. Figure 5 shows the arc branch as it is obtained with values $C=10^{2}$ (bullets, same as before), $C=10^{2}\text{i}$ (crosses), $C=-10^{2}$ (plus signs) and $C=-10^{2}\text{i}$ (squares). All symbols fall onto the same curve, confirming that the position of feedback eigenmodes, which represent singularities in the pseudospectrum, is fixed by the phase of the feedback relation, such that resonance can occur between both ends of the loop.

3 A parallel jet

The effect that imperfect numerical boundary conditions can have on eigenmode computations in a two-dimensional flow domain is now investigated for the case of axisymmetric perturbations in a parallel round jet, governed by the incompressible axisymmetric Navier–Stokes equations,

[TABLE]

Boundary conditions are prescribed as

[TABLE]

The stress-free outflow condition is a common and convenient choice for finite-element computations. Jet radius and centreline velocity are the dimensional length scales in this formulation, and the Reynolds number is chosen as $Re=100$.

The advantage of using a parallel base flow profile is that the global results can be rigorously compared with local instability properties. The standard analytical model of Michalke M84 is adopted,

[TABLE]

and a momentum thickness $\theta=0.1$ is chosen for this example.

The axisymmetric perturbation equations that follow from linearization of (6) about the base flow (8) are discretized with finite elements on a domain of length $x_{max}=20$ and radial extent $r_{max}=50$, using the FreeFEM++ software. This domain is resolved with 80 equidistant elements in $x$, and with 360 non-equidistant elements in $r$. The equations are solved for eigenmodes $[u_{r},u_{x},p]^{T}(r,x,t)=[\hat{u}_{r},\hat{u}_{x},\hat{p}]^{T}(r,x)\exp(-\text{i}\omega t)$, where the eigenvalue $\omega=\omega_{r}+\text{i}\omega_{i}$ contains the angular frequency $\omega_{r}$ and the temporal growth rate $\omega_{i}$.

The global spectrum is shown in Fig. 6. It features a clean upper arc branch (red) with maximum growth rate at $\omega_{1}=1.111-0.052\text{i}$. This mode, labelled ‘1’, is chosen for further analysis. A lower branch (blue) is also present, from which the mode labelled ‘2’ with $\omega_{2}=0.912-1.211\text{i}$ will be examined. Similar lower branches are visible in wake spectra Sipp2007 and in boundary-layer calculations Akervik:2008p1176 .

3.1 Projection onto local instability modes

For comparison, the corresponding local instability problem, for a domain of infinite extent $x\in(-\infty,\infty)$, is solved with a standard Chebyshev collocation technique on a staggered grid Khorrami:1991p1179 , using a coordinate transformation adapted for jet profiles Lesshafft:2007p54 . Identical radial collocation point distributions are used in the global and local computations, such as to eliminate the need for interpolation. Consistent boundary conditions (7c) are imposed in the local problem.

Spatial local instability modes are computed for the $\omega$ values corresponding to the global modes labelled in Fig. 6. The flow is convectively unstable, with an absolute frequency $\omega_{0}=1.074-0.286\text{i}$. Both direct and adjoint modes are solved for. The adjoint local modes represent the dual basis associated with the set of direct modes, and they serve for projecting the spatial structure of the global mode onto the local direct modes.

This projection is carried out in the following way, similar to the procedure used by Rodr\́text{i}guez et al.dani : at a given streamwise station $x$, the radial variations of the global mode perturbation quantities are extracted. Since the eigenvector of the spatial local problem contains auxiliary variables $k\hat{u}_{r}$ and $k\hat{u}_{x}$ (see Lesshafft:2007p54, ), these must also be added to the extracted slice of the global mode. This is accomplished by computing the streamwise derivative of the global $\hat{u}_{r}$ and $\hat{u}_{x}$ fields, and by augmenting the extracted vector $[\hat{u}_{r},\hat{u}_{x},\hat{p}]^{T}$ with $[-\text{i}\partial_{x}\hat{u}_{r},-\text{i}\partial_{x}\hat{u}_{x}]^{T}$. This augmented global slice is finally projected onto the local modes, via scalar multiplication with the associated adjoint modes. As usual, the adjoint modes are normalised beforehand in such a way that their scalar product with the associated direct mode is unity, whereas the direct modes are by themselves scaled to have unit norm.

The spatial local spectrum of the parallel jet profile, shown in Fig. 7 for the frequency values of global modes 1 and 2, is composed of discrete $k^{+}$ and $k^{-}$ modes HM90 , which represent downstream- and upstream-propagating hydrodynamic perturbations inside and near the jet. The labels in Fig. 7 rank all $k^{+}$ and $k^{-}$ modes according to their spatial growth rate $-k_{i}$.

Figure 8 displays absolute values of the projection coefficients, obtained for global mode 1, pertaining to the three dominant local modes. These are the first two $k^{+}$ modes and the first $k^{-}$ mode. Blue symbols represent the $k^{+}_{1}$ mode (the only one displaying unstable spatial growth). Except very near the Dirichlet inlet, the streamwise variation of this local mode amplitude is perfectly exponential, with a spatial growth rate $0.4604$, as measured by a regression fit over $1\leqslant x\leqslant 19$. This value matches within $0.01\%$ the imaginary part of the local eigenvalue. The amplitude of the global mode component $\hat{u}_{x}$ on the jet centerline is shown as a black line for reference.

Green symbols in Fig. 8 denote the amplitude of $k^{+}_{2}$, and red symbols the amplitude of $k^{-}_{1}$. Straight lines indicate the corresponding growth rate of the local eigenvalue for comparison. It is seen that both projections follow the amplitude variations expected from local analysis in a region close to one boundary, where their amplitude is maximal. Clearly, the $k^{+}_{2}$ mode originates at the upstream boundary, whereas the $k^{-}_{1}$ mode is forced at the downstream end. Farther away from those boundaries, both projections approximately follow the slope of the dominant $k^{+}_{1}$ branch. This behaviour results from imperfections in the numerical projection, which apparently only allows a clean distinction between these local modes down to amplitude ratios around $10^{-3}$ in the present setup. Higher spatial resolution in the local and global computations does not improve this threshold. Note that these three modes are highly non-orthogonal, which makes their distinction numerically delicate.

Overall, local mode contributions other than from the $k^{+}_{1}$ mode to the global mode 1 are tractable but rather negligible. Modulations of the centerline velocity perturbation (black line in Fig. 8) are instead attributable to global pressure modes, as will be shown later on.

Very similar results are obtained for the global mode 2, for which the local mode amplitudes are displayed in Fig. 9. However, the local spectrum in this case differs from that of mode 1, as the global frequency $\omega_{2}=0.912-1.211\text{i}$ has an imaginary part below that of the absolute frequency $\omega_{0}=1.074-0.286\text{i}$. Therefore, in the analysis of global mode 2, the local spatial modes are selected from a spectrum where pinching of the $k^{+}_{1}$ and the $k^{-}_{1}$ has already occurred (see Fig. 7b). In this setting, it is now the $k^{+}_{2}$ mode that displays the strongest downstream growth, followed by the $k^{+}_{3}$ mode. A mode from the mixed branch, formed from the $k^{+}_{1}$ and $k^{-}_{1}$ branches after pinching, is denoted $k^{\pm}$.

The global mode 2, as represented by a black line in Fig. 9, is clearly dominated by the $k^{+}_{2}$ mode (blue bullet symbols), but the $k^{+}_{3}$ and $k^{\pm}$ modes are again discernible down to amplitudes three orders of magnitude below $k^{+}_{2}$. It is not clear a priori how the $k^{\pm}$ mode is to be interpreted, in particular with regard to its up- or downstream propagation. However, the projection results plainly show that this mode is generated at the downstream end, from where it propagates upstream; thereby, it behaves as a mode of $k^{-}$ type.

It is stressed that the analysis of spatial branches below the absolute growth rate is indeed meaningful in the present context. Spatial analysis below the absolute growth rate, i.e. after pinching has taken place, is usually said to be in violation of temporal causality, formally expressed by the fact that no integration path can be found in the complex $k$-plane that separates $k^{+}$ and $k^{-}$ branches (Huerre2000, ). This argument however arises in the context of the asymptotic flow behaviour at long times, when indeed the system dynamics are determined by the absolute mode. In the same sense, a global system is asymptotically determined by only the most unstable eigenmode (here: global mode 1). Notwithstanding, global eigenmodes with lesser growth rate do exist, and they are observable in the transient system dynamics. The spatial local modes used in the analysis of global mode 2 are justified in the same way.

In summary, the projections onto local spatial modes demonstrate that global mode 1 is supported by a local $k^{+}_{1}$ mode, whereas global mode 2 relies on a $k^{+}_{2}$ mode. In all likelihood, the same holds true for any global mode of the arc branch ($k^{+}_{1}$) and of the lower branch ($k^{+}_{2}$). Moreover, although not shown in Fig. 6, further lower branches exist at lower growth rates $\omega_{i}<-1.5$, for which a match with higher local $k^{+}$ branches is anticipated. The link between the arc branch and the dominant $k^{+}$ mode has been pointed out in earlier studies (for instance Ehrenstein:2005p1026, ; Akervik:2008p1176, ; Garnaud:2013p1182, ). However, the argument so far has one loose end: the presence of a local $k^{+}$ mode is contingent on it being forced upstream. The present results show that this forcing takes place immediately at the upstream boundary. The essential ingredient that can give rise to a global mode with a $k^{+}$ wave is feedback from downstream.

3.2 Global pressure feedback

The ellipticity of the global linear jet problem is contained in the pressure gradient and in the viscous terms. The latter will only be noticeable over distances much shorter than the numerical box length, and they are not considered in the following analysis. The pressure, as noted by Ehrenstein & Gallaire Ehrenstein:2005p1026 , obeys a Poisson equation, which in the present case of parallel flow takes the form

[TABLE]

with homogenous Dirichlet and Neumann conditions at the upstream and lateral boundaries, respectively, and with $\hat{p}=Re^{-1}\partial_{x}\hat{u}_{x}$ at the outflow.

The following analysis is restricted to the arc branch mode labelled ‘1’ in Fig. 6, but results for mode 2 are not fundamentally different. The pressure amplitude of mode 1 is shown in Fig. 10(a) as $\log_{10}|\hat{p}|$. Its structure is somewhat irregular in the downstream near-field region of the jet, yet the characteristic wavelength $2\pi/k_{r}=4.1$ of the $k^{+}_{1}$ mode is apparent. The stress-free outflow boundary condition is seen to result in $\hat{p}\approx 0$, and the pressure at $r\gtrsim 5$ is essentially a superposition of fundamental solutions of the homogeneous equation $\Delta\hat{p}=0$, with Dirichlet conditions at the inflow and outflow. These are given by

[TABLE]

where $I_{0}$ and $K_{0}$ are the modified Bessel functions of the first and second kind, respectively, and $L=20$ is the streamwise length of the numerical box. The $K_{0}$ functions are exponentially decaying in $r$, manifestly dominant in the present problem, whereas the $I_{0}$ functions grow exponentially in $r$. These only enter the global mode at very low amplitude ($A_{j}/B_{j}\ll 1$) in order to satisfy the Neumann condition at $r_{max}=50$. Only the $\hat{p}_{1}$ component is clearly visible in Fig. 10a, because it experiences the slowest radial decay, but a projection confirms that at least the first five $\hat{p}_{j}$ components enter the pressure field with comparable global norms.

It is known from the analysis in Sect. 3.1 that the global mode involves a strong $k^{+}_{1}$ wave. In the present section, the complementary part is sought that may provide the upstream-reaching part of a feedback loop. Therefore, the $k^{+}_{1}$ component of the global mode is subtracted from the pressure field, using the already known projection coefficients shown in Fig. 8. This ‘stripped’ pressure field $\tilde{p}$ is presented in Fig. 10(b), as $\log_{10}|\tilde{p}|$. A remarkably clean structure is recovered, suggestive of a solution of the Laplace equation $\Delta\tilde{p}=0$ which is forced at the outflow boundary near the jet axis. A small inhomogeneity is also observed at the inflow near the the jet axis.

Note that the local $k^{+}_{1}$ mode by construction represents a particular solution of (9) in a domain of infinite streamwise extent. Within the limits of the simplifying assumption that the propagating perturbations in the region of $\partial_{r}U\neq 0$ are indeed given by the $k^{+}_{1}$ mode alone, in the interior of the bounded domain, the stripped pressure field $\tilde{p}$ only needs to satisfy the inhomogeneities that arise from the boundary conditions. These inhomogeneities on both ends of the domain are thus coupled through the Laplace equation in $\tilde{p}$.

As a final plausibility check, the amplitudes of the supposed downstream-propagating and upstream-reaching components of the feedback loop are compared in Fig. 11 along a path at constant $r$. The radial position of the critical point of the local $k^{+}_{1}$ mode is chosen, $r_{c}=0.92$. A thin blue line represents the pressure amplitude of the $k^{+}_{1}$ mode component, according to the projection carried out in Sect. 3.1, a thick red line represents the amplitude of the pressure feedback signal $\tilde{p}$, and black symbols mark the amplitude of the total pressure field $\hat{p}$. The up- and downstream branches have approximately equal amplitude at the downstream boundary, nearly cancelling each other. Near the upstream boundary, the feedback signal is larger by about a factor 3 than the $k^{+}_{1}$ wave — this appears reasonable in view of the assumption that the latter is forced by the former.

3.3 The influence of box length

Arc branch modes in the literature are consistently found to be sensitive to the length of the numerical box. In some instances (e.g. Garnaud:2013p1182, ), a longer box yields eigenvalues with lower growth rates, in other instances (e.g. Ehrenstein:2005p1026, ) the opposite effect is observed. If the present feedback model is correct, it should allow an estimation of the influence of the numerical box length on the arc branch growth rates.

Eigenmodes of the parallel jet (8) have been computed in numerical domains of streamwise lengths $L=10,15$ and 25, for comparison with the standard configuration $L=20$; in all cases, the radial discretisation and the constant step size $\Delta x$ are unchanged. Resulting modes of the arc branch are displayed in Fig. 12(a). As $L$ is increased, the entire branch is seen to shift to higher growth rates, and the modes are more densely spaced.

The analysis so far suggests that a downstream-convecting $k^{+}$ mode and the elliptic pressure field are coupled in small regions near the inflow and outflow boundaries. In analogy with the Ginzburg–Landau model discussed in Sect. 2, an effective reflection coefficient $C$ may be defined, which relates the forcing of the $k^{+}$ mode at the inflow to the $k^{+}$ amplitude at the outflow, thus lumping the narrow interaction regions into singular actuator and sensor positions in $x$. The effective forcing at the Dirichlet inflow boundary is modelled as

[TABLE]

where $\omega$ is the global eigenmode frequency, and ($k,\hat{q}$) denote the relevant local $k^{+}$ mode. The coefficient $C$ includes the effects of coupling on both ends of the domain, as well as the upstream decay of the pressure signal, seen in Fig. 11.

Numerically obtained values of $C$ are reported in figures 12(b,c), in terms of their absolute value and their phase, for the arc branch in all four numerical domains. In each domain, the reflection coefficient decreases with real frequency, but only weakly so for frequencies larger than one. This effect may be caused by the slower radial decay of $\hat{q}$ at low frequencies, resulting in radially more extended interaction regions on both domain ends. The phase of $C$ varies only slightly with frequency, and it is independent of box length.

If the reflection coefficients are scaled with the cube of the box length, as shown in Fig. 12(d), they neatly collapse onto one curve. Consequently, the effective feedback imparted by boundary reflections decreases with the domain length as $L^{3}$. This scaling factor indicates that the relevant component of the pressure signal, which couples the downstream with the upstream end of the $k^{+}$ wave, is of a quadrupole type (signals from monopole and dipole sources at $x=L$ would decay as $L$ and $L^{2}$, respectively).

It is then straightforward to interpret the effect of an increased numerical box length on a given arc branch mode: if the exponential growth of the dominant $k^{+}$ mode over the added streamwise interval is larger than the algebraic decay of the upstream-reaching pressure signal, over the same added distance, then the spurious inlet forcing will increase in strength, resulting in a higher growth rate $\omega_{i}$. If the $k^{+}$ mode in the added region is locally stable, then the inlet forcing will decrease, and the global growth rate will be lower as a result.

This interpretation appears to be consistent with the two cited examples: the boundary layer investigated by Ehrenstein & Gallaire Ehrenstein:2005p1026 is locally unstable at the outflow, as confirmed by the authors, and longer domain sizes are found to result in higher global growth rates. In contrast, the rapidly spreading jet considered by Garnaud et al.Garnaud:2013p1182 is locally stable at the outflow with respect to axisymmetric perturbations, and indeed lower global growth rates are obtained for the arc branch in longer domains. This criterion may of course be frequency-dependent: certain (complex) frequencies may be locally stable at the outflow while others are unstable, and opposite trends ought to be observed in these frequency ranges. This seems indeed to be the case in the $Re=360$ setting of Coenen et al.Coenen2017 .

4 Possible remedies, tested for a non-parallel jet

While the preceding analyses of the Ginzburg–Landau equation and of the parallel jet served to characterise the global feedback mechanism due to domain truncation, the question how such feedback may be reduced is addressed in this section for the example of a spatially developing jet. With this choice, the results from the previous section can be largely transferred to the new setting.

Each type of open flow may present particular problems in view of domain truncation. In favourable configurations, the numerical boundaries can be placed in stable flow regions, as for instance in a uniform flow upstream of a solid obstacle, like a cylinder. In other cases, the need for truncation may be avoided by prescribing consistent physical boundaries, like confining solid walls, which are then part of the flow configuration.

Jets belong to a more problematic category. They are created by some upstream source of momentum, and it is in general not feasible to include the entire upstream apparatus (fan, chamber, nozzle, etc.) in the calculations, which would furthermore defeat the purpose of any generic description of jet dynamics. The formulation of upstream flow and boundary conditions is therefore necessarily imperfect with respect to any flow realisation, as potentially important regions are not accounted for.

For the present study, the jet is modelled as issuing from an orifice in a solid wall, very similar to the configuration of Garnaud et al.Garnaud:2013p1182 . Upstream of the orifice, the flow develops in a straight circular pipe, from where it exits with a fairly thin boundary layer. The base flow is computed as a steady, axisymmetric solution of the incompressible Navier–Stokes equations (6), with an inflow condition

[TABLE]

imposed at $x=-10$. Newton–Raphson iterations are performed in order to converge to a steady flow state. The numerical domain for these base flow calculations is truncated downstream at $x_{max}=50$ and radially at $r_{max}=50$, where stress-free conditions are applied. The Reynolds number is set to 1000 in the free jet. However, in order to maintain a thin shear layer at the orifice, this value is varied exponentially inside the pipe, between $Re=10^{5}$ at $x=-10$ and $Re=10^{3}$ at $x\geqslant-0.2$, in the base flow calculation. The resulting flow field near the orifice at $x=0$ is represented in Fig. 13, together with the characteristic streamwise variations of the shear layer momentum thickness $\theta_{m}$ and of the centreline velocity.

Eigenmode calculations are carried out on smaller domains, where portions of the base flow are cropped at the inflow, the outflow, or both, in order to probe the effect of domain truncation. All configurations are listed in table 1. The Reynolds number in these calculations is set to 1000 throughout the domain.

The influence of upstream truncation on eigenmodes is considered first. Figure 14a shows spectra for three different domains, with homogeneous Dirichlet conditions (7a) imposed at $x_{min}=(0,-5,-10)$, respectively. The downstream end of the numerical domain is placed at $x_{max}=40$, where stress-free conditions (7b) are applied in all three cases. Unconverged and lower-branch eigenvalues are not shown here and in the following for clarity; the criterion for convergence is that an eigenvalue could be reproduced within three-digit accuracy using two different shift values.

It is seen from Fig. 14a that the arc branch is quite insensitive to the position of the upstream Dirichlet boundary. Even a full truncation of the upstream pipe (case 1) only results in a slight stabilisation at high frequencies. Eigenvalues obtained with pipe lengths of 5 and 10 radii (cases 2 and 3) are virtually identical. Stress-free inflow conditions (cases 4, 5 and 6, spectra shown in Fig. 14b) are found to perform less favourably. When applied at $x=0$, these conditions give rise to global instability; the inclusion of portions of the pipe has a stabilising effect, but even with $x_{min}=-10$ the growth rates along the entire arc branch are still significantly higher than those obtained with Dirichlet conditions.

It is to be expected that global pressure fluctuations induce a generation of vorticity waves at the solid corner at $x=0$, and the independence of the eigenvalues in cases 2 and 3 (Fig. 14a) on the upstream boundary position suggests that this physical effect is dominant over spurious pressure-vorticity coupling at $x_{min}$. It remains to be determined to what extent the downstream boundary conditions emit spurious pressure signals, and how these may be reduced.

Garnaud et al. Garnaud:2013p1182 tested stress-free and convective outflow conditions in an almost identical jet configuration, and found no significant difference in the resulting spectra. In the present study, the potential of absorbing layers (C04, ) for a stabilisation of the arc branch is examined. To this end, an artificial damping term is added to the linear perturbation equations, with the purpose of reducing perturbation amplitudes before they reach the numerical boundary at $x_{max}$. If the original global eigenvalue problem is written

[TABLE]

then the eigenvalue problem with absorbing layer is defined as

[TABLE]

For cases 7 and 8 (see table 1), the damping parameter is prescribed as

[TABLE]

This variation of $\lambda$ is continuous in its first derivative, and it has maximum values 0.5 and 0.75 at $x_{max}=40$ and $x_{max}=50$, respectively.

The effect of absorbing layers on the eigenvalue spectrum is documented in Fig. 15. Damping in the downstream region $19\leqslant x\leqslant 40$ (case 7) is seen to reduce the growth rates of arc branch modes significantly (compared to case 2). Eigenfunctions of selected modes, marked by circles in Fig. 15, are compared in Fig. 16. An integral amplitude measure is defined as

[TABLE]

and the amplitude curves in Fig. 16 are normalised with respect to their values at the outflow boundary for comparison. While the red curve (case 2, no absorbing layer) shows a monotonic growth of perturbation amplitude throughout the domain, similar to the parallel jet results in figures 8 and 9, the blue curve (case 7) reaches a maximum inside the absorbing layer and subsequently decays in $x$. Remarkably, the ratio $A(x_{max})/A(0)$ is identical in both cases, which is consistent with the interpretation of these modes as being the result of spurious feedback from the outflow, in the same way as detailed in Sect. 3 for the parallel jet. The absorbing layer reduces the gain of the $k^{+}$ hydrodynamic branch of the feedback loop, and thereby the temporal modal growth rate.

Based on the discussion in Sect. 3, stronger artificial damping of the hydrodynamic branch ought to lead to ever smaller modal growth rates; longer domains should furthermore lead to a reduced amplitude of the incident pressure signal upstream, due to its cubic decay. Both conditions are combined in case 8 (see table 1), of which results are included in figures 15 and 16. Indeed, the temporal growth rates are further reduced, and the spatial decay of the eigenfunction amplitude in the absorbing layer is more pronounced than in the case 7.

However, the shape of the arc branch in case 8 displays some differences with respect to all other cases, and similar mode patterns have been reported from jet calculations on long domains by Coenen et al. Coenen2017 , where no artificial damping was applied. In the vicinity of the least stable mode (circled in Fig. 15), a regular spacing in $\omega_{r}$ is still observed, but with larger distances between consecutive modes than in cases 2 and 7. The corresponding amplitude function in Fig. 16 shows that the ratio $A(x_{max})/A(0)$ is smaller than in the two other cases, whereas the cubic decay of the reflection coefficient with box length, as described in Sect. 3, should have resulted in a larger ratio. Furthermore, stronger damping through higher values of $\lambda$, tried in test calculations that are not shown here, does not decrease the maximum growth rate much further, but it quickly leads to ill-conditioned system matrices. Case 8 seems indeed to mark an efficiency limit of absorbing layers for the stabilisation of the arc branch.

The above observations, in particular the increased spacing of modes in $\omega_{r}$, suggest that pressure feedback in the strongly damped case 8 may originate from the interior of the domain, possibly from the location of the amplitude maximum. Such feedback may be spurious, in the sense of Heaton et al. heaton , or it may be physical. This hypothesis is suggested for further examination.

5 Conclusions

A family of linear instability eigenmodes, named the arc branch, has been analysed in view of its physical or numerical origin. Branches of this type have previously been observed in a large variety of open flows. All results presented herein lead to the conclusion that arc branch modes in incompressible flow calculations are an artifact of domain truncation, due to spurious pressure feedback between numerical inflow and outflow boundaries.

A Ginzburg–Landau model was investigated first, because in this setting the effect of explicitly prescribed feedback between a downstream sensor and an upstream actuator could be examined without ambiguity. In the presence of feedback, a branch of eigenmodes was observed to arise that exhibits all typical characteristics of the arc branch described by Coenen et al.Coenen2017 for the example of light jets. These modes align with a contour of the pseudospectrum, spaced at regular intervals of the real frequency, and their spatial eigenfunction is characterised by an integer number of wavelengths between actuator and sensor locations. As the strength of the feedback is increased, the branch moves steadily upward to higher growth rates in the complex frequency plane. Eigenmodes of the feedback-free system that lie below the arc branch are no longer detectable. An important observation is that this arc branch in the Ginzburg–Landau equation with feedback has no counterpart in the feedback-free spectrum: these eigenmodes are not merely affected by the presence of feedback, indeed without it they do not exist.

An incompressible parallel jet in a truncated domain was examined next. No explicit feedback was prescribed in this setting, but a similar arc branch of eigenmodes was nonetheless found to dominate the spectrum. A lower branch of eigenmodes with stronger temporal decay was also described, which may correspond to subdominant branches observed in wakes and boundary layers. Both the arc-branch and the lower-branch modes in the parallel jet were shown to be composed primarily of one downstream-propagating $k^{+}$ wave, as computed from a local spatial analysis. Each of the two branches involves a different $k^{+}$ mode. This observation raises the question how a $k^{+}$ wave is generated at the upstream boundary of the numerical domain. The analogy with the Ginzburg–Landau model from Sect. 2 suggests feedback from downstream.

The most plausible mechanism for global feedback in a truncated domain is the generation of pressure perturbations due to the artificial downstream boundary condition, as described by Buell & HuerreBuellHuerre . The pressure field of the least stable arc branch mode was decomposed into one portion associated with the prominent $k^{+}$ wave, which cannot be involved in upstream feedback, and into a residual part that is essentially governed by a Laplace equation forced at the domain boundaries. The latter appears to be strongly dominated by a source region at the outflow, near the jet axis. A numerical evaluation of effective reflection coefficients, in analogy with the Ginzburg–Landau model of Sect. 2, established as a principal result that the strength of upstream feedback decays with the numerical box length to the third power. This scaling is indicative of a spurious pressure quadrupole situated at the outflow. The algebraic nature of the pressure decay allowed a prediction of the effect of box length variations on arc-branch growth rates, depending on the local stability or instability of the flow, apparently consistent with common observations. While this study has been limited to incompressible flow settings, similar mechanisms will also be present in compressible calculations, with the difference that spurious pressure signals are then propagated by a wave equation. This may in fact result in stronger feedback, as the far-field acoustic pressure only decays with the first power of the distance from the source.

Finally, possible strategies for a reduction of spurious pressure feedback have been examined for the example of a spreading jet. On the one hand, upstream boundary conditions must be chosen that minimise the unphysical conversion of pressure feedback from downstream to vortical perturbations. It has been found, in this particular flow example, that Dirichlet conditions perform much better than stress-free conditions in this regard. On the other hand, the generation of spurious pressure signals from the outflow boundary must be reduced. It has been demonstrated that artificial damping in a downstream absorbing layer provides an efficient means to achieve this. An advantage of this technique is that it is straightforward to implement in any usual open shear flow problem. It has been noted, however, that increasingly strong damping does not reduce the growth rate of the arc branch to arbitrarily low levels.

Acknowledgements.

Discussions with Xavier Garnaud and Wilfried Coenen greatly helped shape the ideas presented in this paper. The study was supported by the Agence Nationale de la Recherche under the Cool Jazz project, grant number ANR-12-BS09-0024.