Qualitative dynamics of wavepackets in turbulent jets

Qualitative dynamics of wavepackets in turbulent jets

Onofrio Semeraro11footnotemark: 1 semeraro [at] ladhyx.polytechnique.fr François Lusseyran22footnotemark: 2 Luc Pastur33footnotemark: 3 Peter Jordan44footnotemark: 4 LadHyX, CNRS, École polytechnique – Palaiseau, FR LIMSI – CNRS, Orsay, France Institut PPRIME, CNRS, Université de Poitiers, ENSMA – Poitiers, FR

It has long been established that turbulent jets comprise large-scale coherent structures, now more commonly referred to as “wavepackets” (jordan2013wave, ). These structures exhibit a remarkable spatio-temporal organisation, despite turbulence.

In this work we analyse, from a qualitative point of view, the temporal dynamics of axisymmetric wavepackets educed, experimentally, from subsonic iso-thermal jets. We use the data presented by breakey2013near (), where time-series of the wavepackets are extracted at different streamwise locations. A thorough analysis is performed; statistical tools are used for estimating the embedding and correlation dimensions characterising the dynamical system. System identification is used for computing nonlinear surrogate models. Finally, control-oriented linear models are computed.

The goal of the contribution is to assess the extent to which non-linear models are necessary, or appropriate, for description of the temporal wave-packet dynamics and to provide a complementary perspective to the current modelling.

journal: Physical Review Fluids (PRF)

1 Introduction

Following the early studies of Mollö-Christensen (mollo1960sound, ; mollo1964experiments, ; mollo1967jet, ), a considerable body of work has been devoted to exploring the nature of organised motions that are observed in turbulent jets (Crow-Champagne, ; lau1972intrinsic, ; Michalke-JFM-1975, ; armstrong1977coherent, ; moore-JFM-1977, ). It was hypothesised early on that this component of the flow might be understood in terms of an instability of the turbulent mean (Crow-Champagne, ; michalke1971instability, ; bechert1975amplification, ; cohen1987evolution, ; crighton1976stability, ; plaschko1979helical, ); and the importance of such flow structures for sound radiation has been suggested in numerous studies (michalke1970wave, ; michalke1971expansion, ; bishop1971noise, ; Crow-PhysSoc-1972, ; Crighton-Theory-1975, ; Michalke-JFM-1975, ; crighton1990shear, ).

It is only more recently, however, thanks largely to progress in theory, numerical simulation and experimental diagnostics, that it has been possible to explore these hypotheses in a comprehensive manner. Exhaustive comparison of the results of theory with experimental measurements has confirmed that the average characteristics of coherent structures in turbulent jets are remarkably well described by solutions of the Navier-Stokes equations linearised about the turbulent mean (suzuki2006instability, ; cavalieri2013wavepackets, ; jordan2014modeling, ); these solutions are synonymous with globally stable modal solutions (nichols2011global, ; garnaud2013preferred, ; schmidt2016super, ), which, physically, amount to hydrodynamic waves that are convectively amplified in the upstream region of the flow, but become neutrally stable and then decay farther downstream. It is for this reason that they have been given the denomination ‘wavepacket’ (jordan2013wave, ). An example of a modal solution, from schmidt2016super (), is shown in figure 1, where it is compared with a realisation taken directly from the Large Eddy Simulation (bres2015large, ) that provided the mean flow. The wavepackets can be seen to comprise a hydrodynamic component, with an amplitude envelope as described above, and an acoustic component that takes the form of a directive beam radiating to shallow angles.

Figure 1: Left: realisation taken from Large Eddy Simulation (axisymmetric component of pressure at ). Right: Global mode . Courtesy of schmidt2016super ()

Considering the Reynolds number of the jet, , and the fact that it issues from a nozzle with fully turbulent boundary layers, the agreement between the LES realisation and the linear global mode is striking. One is compelled to ask how, despite the non-linear, orderless character of the turbulence that dominates the fluctuation energy of this flow, such organisation can exist and how a linear model can capture so many of its features.

However, something less clear from the images is the considerably lower acoustic efficiency of the linear wavepacket. As shown in a number of recent studies (zhang2014just, ; jordan2014modeling, ; baqui2015coherence, ), the intensity of the soundfield of the linear solution, once the hydrodynamic fluctuation levels of LES and model have been matched, can be as much as 40dB below that of the LES. Such a large discrepancy, in the face of such compelling organisational similarity, is intriguing and the subject of a number of ongoing research efforts, which we summarise briefly in what follows.

The present understanding is that the linear dynamics do not contain the acoustically important degrees of freedom: the average wavepacket does not generate the average sound (cavalieri2011jittering, ; cavalieri2011using, ; kerherve2012educing, ; jordan2014modeling, ; zhang2014just, ). The essential sound-producing motions are those associated with higher-order statistics of the wavepacket dynamics. These motions have been denoted ‘jitter’ by (cavalieri2011jittering, ), and would appear to be underpinned by non-linearity (jordan2014modeling, ; zhang2014just, ; towne2015stochastic, ; tissot2016sensitivity, ).

The nature of this non-linearity remains to be clarified, and it motivates the present study. Do wavepackets jitter on account of non-linear wave-wave interactions (sandham2006nonlinear, ; suponitsky2010linear, ), i.e. is the non-linearity of an intrinsic kind? Or is it rather an extrinsic, stochastic forcing of linear waves by background turbulence that leads to the activation of these higher-order degrees of freedom (towne2015stochastic, ; tissot2016sensitivity, )? Or might both extrinsic and intrinsic mechanisms be at work? In this work, in light of the foregoing considerations and questions, we are interested in low-dimensional projections of the complete wavepacket dynamics.

A further issue, that motivates both the modelling efforts evoked above and the work we undertake in this paper, is that of control. Given the success of linear models in predicting both the average wavepacket structure, in frequency space, and the real-time evolution in an experimental context sasaki2015real (), is it legitimate to ask if the tools provided by linear control theory might be sufficient? Such possibilities are presently being explored by sasaki2016closed (). But if it were to be necessary to directly manipulate the more subtle dynamics associated with the energetically-unimportant-but-acoustically-essential degrees-of-freedom, then a non-linear control framework would be required (machine learning etc. gautier2015closed ()). In both cases the question of the dimension of the space spanned by the dynamics must again be asked: is this small enough for control to be realistic in an experimental context?

And, finally, there is the issue of the physical interpretation of POD modes, which some recent results show may correspond to the highest gain structures of the linear Navier-Stokes operator subject to stochastic external forcing by turbulence. With this in mind, the various analyses of the dynamics will be performed using both raw and POD-filtered data.

Consideration of the above issues via qualitative dynamical analysis has, to the best of our knowledge, not previously been attempted. We therefore consider the same turbulent jets studied by (cavalieri2012axisymmetric, ; cavalieri2013wavepackets, ; breakey2013near, ), using a variety of non-linear signal-processing and system identification techniques, in order to explore dynamics of low-frequency wavepackets. Particular attention is given to pressure signals dominated by fluctuations at Strouhal numbers , where linear models are found to perform most poorly.

Our analysis starts with a minimal description of the experimental setup in Sec. 2; we refer to breakey2013near () for a detailed description. The remainder of the paper is organised in two parts. The behaviour of the time-series obtained from the experimental campaign is discussed in Sec. 3, where statistical tools and spectral analysis are introduced for a preliminary characterisation of the temporal dynamics.

In the second part, we pursue an alternative approach by analysing reduced-order models of the time-series, based on system identification; the approach is introduced in Sec. 4, while the dynamical analysis is reported in Sec. 5 and Sec. 6 . More details on the techniques are given in A. The manuscript finalizes with a discussion of the results and concluding remarks in Sec. 7.

2 Experimental data

The study presented in this paper is based on the analysis of the temporal behaviour of axisymmetric wavepackets educed from iso-thermal turbulent jets. In particular, we consider the near-field pressure fluctuations measured by breakey2013near (). For a detailed description of the measurement setup, data acquisition and post-processing the reader can refer to that paper. We here limit ourselves to a brief presentation of the data, in order to highlight those features, which are relevant to the analysis we perform.

Figure 2: Overview of the experiment, adapted from breakey2013near ().

2.1 Experimental setup

We consider an unforced, isothermal, subsonic jet issuing from a round nozzle with a fully turbulent boundary layers. The flow is thus very different to the transitional jet studied by narayanan2013dynamical (). The setup is shown in Fig. 2, and a schematic of the near-field microphone array used for eduction of the axisymmetric wavepacket signature is sketched in Fig. 2. The experiments were carried out in the anechoïc free jet facility at the Centre d’Etudes Aérodynamiques et Thermiques (CEAT), Institut Pprime, Poitiers (FR). The jet Mach number ranged from to . In what follows, the diameter of the nozzle, , is taken as the reference length-scale, while the exit velocity, , is the reference speed. Time is scaled as . The corresponding Reynolds number ranges between , where is the density and the air viscosity. The potential core length of the jet ranges between .

Figure 3: The near-field signal is shown as a function of the distance from the nozzle exit () and the time , scaled with the velocity and the diameter . The measured quantity is the pressure related to the axisymmetric mode . The time-series are obtained from the POD reconstruction described in breakey2013near ().

Azimuthal rings of six microphones recorded pressure fluctuations from the near-field of the jet, on a conical surface. An azimuthal ring of three microphones recorded the far-field fluctuations. In this work, we focus on the near-field fluctuations; two different datasets are analysed, obtained from two measurement campaigns.

  1. The first set of data is obtained by simultaneous measurements at different streamwise locations in the near-field. At each location, an azimuthal ring, each with microphones, is positioned. The array thus comprises a total of microphones, placed on a conical surface in the near-field of the jet, as shown in Fig. 2. The spacing between the rings is . At each location a Fourier-series decomposition is performed in the azimuthal direction.

  2. The measurements of the second campaign are performed using a rig of equi-spaced, azimuthal rings, which are displaced in the streamwise direction over the range with a streamwise resolution of . This set of measurements allows a more detailed computation of the two-point flow statistics, providing the cross-spectral-density matrix, an eigen-decomposition of which provides POD modes.

This POD basis is used for the projection of the data from the first campaign. In this way, a larger domain along is spanned. The temporal dynamics of the axisymmetric mode is shown in Fig. 3, as a function of time and the streamwise direction. This space-time picture is complementary to the frequency-space realisation shown in Fig. 1. Both show clearly the highly organised nature of wavepackets in these high-Reynolds-number, fully turbulent jets.

Our study focuses on the temporal dynamics of the axisymmetric mode . The analysis is carried out by considering the time-series extracted at given streamwise locations, for each of the datasets. By definition, a time-series is a sequence of equi-spaced data points, function of time. For our case, the sampling rate is , and the total number of points for each series is .

In the following section, preliminary analysis of the temporal behaviour is performed by means of spectral analysis and estimation – using non-linear tools – of the embedding and correlation dimensions associated with these wavepackets.

3 Part I: Time-series analysis

In this section, we consider both the original data and the time-series post-processed by means of proper-orthogonal decomposition (POD). In particular, we study the runs at and Reynolds number , at different locations along the streamwise direction (, ).

The main goal in this section is to characterise the time–series from a dynamical-system point of view. Although we are considering a fully turbulent case, we study a projection of the system on a set of linear, axisymmetric wavepackets. Thus, the working hypothesis is that the resulting temporal dynamics is low-dimensional. In principle, a dynamical system can be qualitatively analysed by reconstructing the – unknown – underlying attractor from a sequence of observables of its state, for instance the time-series. This idea is at the heart of the Takens’ theorem takens1981detecting (), a delay embedding theorem providing the conditions under which the reconstructed dynamics preserves the properties of the original system. In particular, the theorem ensures a diffeomorphism between the original – unknown – phase space and the reconstructed phase space. In that sense, there are two crucial parameters: i) the minimal embedding dimension , i.e the dimension of the projection phase space of the reconstructed dynamics; ii) the correlation dimension , a measure of the fractal dimension of an attractor (see grassberger1983characterization ()). According to these values one can evaluate the dimensions of the attractor of the system and – if possible – study its projection.

In general, the available time series are univariate; thus, the embedding is typically the first step for the qualitative analysis of a dynamical system. Different procedures of embedding are available. For instance, a singular-value decomposition of the Hankel matrix based on univariate time-series can be performed, where the dimensions of the rectangular matrix are the length of the time series and the maximum embedding. A second alternative consists of computing higher-order derivatives of the time-series, such that each of the derivatives corresponds to one of the columns of the embedding vector. In this work, we rely on the most classical approach: once a maximal embedding dimension is chosen, each coordinate is obtained by a time shift ; thus, introducing the time series , the embedding vector is obtained as


where indicates the time-span; note that the last channel will have a total shift of with respect to the first one. Once an embedding vector is obtained, we need to identify the minimal embedding dimension that allows determination of the dimension of the (hypothesised) attractor. A small indicates the possibility of unfolding the attractor. A popular method for estimating the minimum embedding dimension is the false nearest algorithm, kantz2004schreiber (); cao1997practical ().

Figure 4: Temporal behaviour for the mode at . In the inset the original signal (TS-Raw) is compared to the result of the POD post-processing (TS-POD), at a location . A second non-linear filter is applied to both the series as shown in the insets and , where the original time-series (black solid line) is compared to the filtered one (red dashed line); the resulting error is shown as a blue dashed-dotted line. Finally, in the inset the spectrum of the four time-series is shown. A difference between TS-Raw and the TS-POD is observed for , while good agreement is obtained at numbers relevant for our analysis.

The second dimension of interest is the correlation dimension , providing an estimate of the fractal dimension of an attractor. This quantity is related to the embedding dimension by the Takens’ theorem, stating that (takens1981detecting ()).

In what follows, we adopt statistical tools described from the theoretical point of view in the book by Kantz & Schreiber kantz2004schreiber (). A thorough discussion of these techniques is beyond the scope of this paper: we will focus mainly on the validation of the results. Unfortunately, these statistical methods are capable of producing results also in the presence of purely random time-series. Thus, the validation has a twofold goal: assessing the robustness of the results, and confirming the determinism of the dataset.

3.1 Filtering and spectral analysis

Figure 5: Spectral analysis for the series based on POD projection, at different locations along , namely , , and . It can be observed that the peak shifts from to when considering downstream locations.

In order to analyse the time-series, we first perform non-linear noise reduction. Indeed, noise is a dominant, limiting factor for the embedding procedure kantz2004schreiber (). The de-noising algorithm is a moving-average filter, applied along the trajectory identified in the embedding space, on the assumption that the dynamics is continuous. In Fig. 4, the time-series is shown at . In particular, we are interested in assessing to what extent the applied filters modify the essential dynamics of the time-series. In the following, for ease of discussion, the original time-series is denoted TS-Raw while data filtered using POD as the projection basis is denoted TS-POD.

In Fig. 4, the TS-Raw data are compared with the TS-POD set in the time-domain. A quantitative assessment is reported in Fig. 4 using Welch’s power spectral density (PSD) estimate. At the significant Strouhal numbers, , it can be see that the POD filter does not strongly modify the dynamics of the time-series. Differences can be observed only at . At lower , we observe the cut-off due to high-pass filtering of the original data applied to remove energy below the anechoïc cut-off frequency of the windtunnel. The two data sets are filtered using the non-linear filters in Fig. 4 and Fig. 4. By definition, the non-linear filter does not act on specific bandwidths, but on the whole spectrum, while preserving the foliation of the attractor in the embedded space. Welch analysis confirms this behaviour: in Fig. 4 we can observe that the non-linear filter acts on the whole range of numbers. The standard deviation of the removed noise is and , for TS-Raw and TS-POD, respectively.

Figure 6: Spectrogram for the time-series based on the POD projection at . In the insets and , the power spectrum is shown as a function of the Strouhal number at and .

The spectral content of the TS-POD data at is compared in Fig. 5. We observe that the maximum fluctuation levels move progressively to lower frequencies as the observation position moves downstream, from at to at . A change in the maximum amplitude can also be observed. An alternative analysis is provided by the spectrogram, shown in Fig. 6 for the TS-POD at . The graph is in two dimensions: the horizontal axis represents the time scale, while the number is reported along the vertical axis. The amplitudes are shown as a colour map: at each instant, the intensity of the fluctuation is shown as a function of . The Welch estimate can be regarded as an average along the time-span of the results shown in the spectrogram. In particular, in Fig. 6, two instants are shown at and , respectively. The instantaneous amplitudes of certain frequencies in the bandwidth can be observed; this behaviour is a time-frequency view of the organised behaviour manifest in space-time in Fig. 3. These behaviours are the signature of the amplifier nature of the flow: wavepackets are hydrodynamic instability waves that exist on the turbulent mean flow. These waves acquire initial amplitudes and phase at upstream stations – either from disturbances issuing from the turbulent motions within the nozzle, or from distributed turbulence that acts as a volume force – and evolve in the downstream direction according to the stability properties of the linear operator, possibly subject also to distributed forcing from background turbulence; this last issue is one that is presently being studied.

3.2 Embedding dimension and correlation dimension

Assessment of the correlation dimension, , and embedding dimension, , is performed for both TS-Raw and TS-POD datasets, filtered non-linearly. The results are collected in Fig. 8, as a function of the streamwise direction .

3.2.1 Minimal embedding dimension

The minimal embedding dimension is estimated by using the algorithm described by Cao in cao1997practical (). The algorithm does not strongly depend on the length of the time series and it is capable of estimating the embedding dimension also for time series describing high-dimensional attractors. More importantly, it allows to clearly determine whether the signal is deterministic or stochastic, as shown in Fig. 7. The dimension is estimated by analysing the behaviour of the parameter : the saturation indicates the minimal embedding dimension. We chose it using a threshold value , giving . The quantity explicitly accounts for the correlation: random series will be characterised by , for any . For deterministic data, cannot be a constant. In our case, we find that the values of are not constantly unitary, as shown in Fig. 7: this is a first clue that we are analysing deterministic time series.

Figure 7: Example of minimum embedding dimension estimation. The POD-filtered data are used, in . The saturation of the parameter (blue curve) indicates the minimal embedding ; the complementary parameter (red curve) confirms the determinism of the dataset, as for .

A critical aspect is the choice of the embedding delay . In principle, the embedding dimension is independent of the time delay; in practice, the minimum embedding dimension estimate may depend on this choice and needs to be verified case by case. In Fig. 8, we identify a confidence interval for the estimates; the upper bound corresponds to a time-delay , while the lower bound is obtained for . The choice of the embedding delay is done by analysing the auto-correlations. points have been used for the computations, although we verified that is scarcely affected by this choice.

The confidence intervals are shown using shadowed areas. Each point corresponds to the value where the parameter saturates, as discussed in Fig. 7. Dotted curves identify the minimal embedding dimension computed by averaging the results obtained with points and . It can be observed that the TS-Raw series is characterised by a rather constant value of , oscillating between a minimum value at and maximum values further upstream, at . This estimate follows closely the upper bound of our confidence interval. The lowest value is found at , for the runs at .

The behaviour of the TS-POD series is coherent with the original data, although we found higher values. This apparently counter-intuitive behaviour is discussed later (see Sec. 3.2.3). For values , the minimal embedding dimension is practically constant until ; however, note that in this region the data are extrapolated, so relevant dynamics carried by higher order POD modes is missing.

3.2.2 Correlation dimension

Figure 8: The estimates of the embedding dimension and the correlation dimension are summarised as function of for both the datasets under investigation. The minimal embedding dimension is computed for two different time-delays, defining the confidence intervals indicated with shadowed areas: a smaller time of embedding, , defines the upper limit, while a longer time, , is related to the lower bound. The red area is obtained for the TS-Raw series, while the grey one is associated with the results of the dataset TS-POD. The nominal estimate for (dotted curves) is obtained by averaging the results of different runs. The correlation dimension is indicated with a squared-dotted black curve and a squared-dotted red curve for TS-POD and TS-Raw, respectively. A green line indicates the Eckmann-Ruelle limit: values of above this bound are not physical.

The correlation dimension is estimated using the routines included in hegger1999practical (), based on the Grassberger-Procaccia algorithm kantz2004schreiber (). Convergence tests were performed over the input parameters of the algorithm: the embedding quantities , the time delay, and the Theiler window. The Theiler window accounts for the possible oversampling. In fact, the computation of is based on a correlation sum: in the time-series, successive elements are generally not independent and can be highly correlated. To limit this effect, that could lead to inconsistent results, a time-shift – i.e. the Theiler window – is introduced to reduce the correlation between points during the pair counting.

A limit characterises the computation of the correlation dimension; in particular, the value of the correlation dimension over a decade cannot exceed , where is the number of points in the time series (Eckmann-Ruelle limit, eckmann1992fundamental ()). The amount of data available allows us to get a value of . This means that when we approach this upper bound, it is not legitimate to conclude that the value corresponds to the dimension of the inner dynamics.

In Fig. 8, the results are shown with squared-dotted curves. For each positions of the available data-series, we analysed two subsequent blocks of data containing a total of points; thus, only values of correlation dimension are considered relevant (a green shadowed area indicates the limit). The final value of is obtained by averaging the results of the blocks. We observe that the time-delay imposed for the embedding procedure does not influence the final result. A Theiler windows of was chosen; we do not observe relevant changes when increasing this parameter.

For the series TS-Raw, we observe a behaviour similar to , with in the span . For the TS-POD, consistent with the previous estimate of , we found that the value of are slightly higher. In the upstream region , the value close to the upper bound does not allow us to conclude that this is the effective value of the correlation dimension. However, the dimension with respect to the embedding dimension is bounded within the limits imposed by the Takens theorem.

A decrease is observed at the location extrapolated in the vicinity of the nozzle and far downstream, for locations . Also in this case, low values observed in the TS-POD data are due to the POD filtering suppressing important dynamics in these regions of the flow (see breakey2013near ()).

3.2.3 A brief discussion on the estimated dimensions

(a) Phase space,
(b) Phase space,
Figure 9: Phase-space analysis for two positions taken from the TS-POD series. The embedding space is obtained by time-delayed coordinates, using .

The preliminary analysis of the minimal embedding dimension and the correlation dimension indicates deterministic behaviour, despite the limitations that both algorithms pose. Interestingly, the values of and are higher when considering the POD-filtered data. This result has been found also by computing the embedding dimension with a different algorithm based on the false nearest strategy (results are not reported). We believe that the behaviour is due to the physics of the system described in the different datasets and it is not a numerical artefact. As mentioned in Sec. 2, the simultaneous pressure measurements of TS-Raw are projected on a basis of PODs obtained using two-point flow statistics on a larger number of points in the streamwise direction leading to a finer resolution. Although from the spectral point of view the resulting time-series are equivalent, except at higher (see Fig. 4), from the dynamical point of view the projection on the POD-basis might lead to a richer dynamics due to the increased accuracy of the measurements. When considering the dimensions outside the range , i.e. where the simultaneous measurements are not available, we observe a rather low correlation dimension, while the embedding dimension is constant until . This effect is due to the extrapolation of the data in this region of the flow. For this reason, in the following sections we focus only on the points .

(a) Phase space,
(b) First return map,
Figure 10: Phase-space analysis for the TS-POD series at (left) and first return map obtained by stacking the local peaks of the time series. The embedding space is obtained by time-delayed coordinates, using .

The embedding space that should be used according to the Takens’ criterion is rather large. In Fig. 9 and Fig. 10, we show three phase-portraits reconstructed in a three dimensional embedding space obtained by imposing a delay of , between each of the three coordinates. In all cases, at a first glance, the reconstructed phase space suggests a toroidal nature for the dynamical system; unfortunately, the projection is poor due to the high-dimensionality and cannot be further investigated. Moreover, the first return map in Fig.  suggests contamination of the dynamics due to external forcings. Similar results were obtained for all the positions along regardless of which embedding strategy was used (SVD embedding and derivative embedding).

In conclusion, the results are consistent with a determinism of the time-series under investigation and suggest that the dynamics is relatively low-dimensional. However, the embedding and correlation dimensions are still too large for useful characterisation of the attractor geometry. For this reason, alternative strategies are necessary for exploration of the dynamics of the system. In the next section, system identification is applied: with the same tools, we attempt to replace the original data with surrogate models that capture a part of the dynamics.

4 Part II: A system-identification modelling approach

In the previous section, the preliminary study of the time-series – based on statistical tools and Fourier analysis – provided clues that the underlying temporal dynamics of the axisymmetric coherent structures is deterministic. However, the dynamics educed locally is rather rich and cannot be analysed by simply using the standard tools from dynamical system theory. From a global view point, the turbulent jet is characterised by a strong convective behaviour: the intrinsic dynamics is sensitive to external forcing huerre1990local (). Thus, one might wonder to what extent the behaviour educed from the time-series analysis is due to intrinsic non-linear interactions, external activation of degrees of freedom driven by external (stochastic) forcings or a combination of intrinsic and extrinsic mechanisms.

In an effort to understand the underlying dynamics of the system, surrogate models can be introduced by applying system identification kantz2004schreiber (); aguirre2009modeling (). By definition, system identification aims at building mathematical representations of dynamical systems from set of observables, i.e. measured data. In other words, given a dynamical system, we can analyse it from its observables and infer a model reproducing the dynamics (or part of it). Starting from the hypothesis that the chosen observables properly represent the system behaviour, this approach can be summarised in the following steps:

  1. Identification of the model structure and parameter estimation – models are identified such that the dynamics of the observables is reproduced. Due to the number of parameters involved, the number of possible models is typically large (see Sec .4.1 and A).

  2. Model validation – among all the possible models, only a few allow a proper characterisation of the dynamical system. Validation is necessary in order to discard the models unable to reproduce a dynamics close to the real one.

  3. Prediction or dynamical analysis – according to the goal of the identification process, we can aim at minimising the prediction error given a time horizon or analysing the dynamical properties of the system.

The last aspect is particularly relevant: according to the goal of the identification process, the choices of the inputs/outputs as well as the model structure are of great importance.

In the following, we will often refer to noise indicating with this term the external forcings driving the system dynamics. As already said, turbulent jets are convectively unstable, so strongly sensitive to external forcings; measurement noise will be generally indicated as error or residual in the context of the identification problem.

4.1 Choice of the model, inputs and outputs

For our application, the chosen model consists of non-linear Volterra series of polynomials. Hereafter, the outputs will be denoted by and will be represented by the measurements taken at different locations along the streamwise direction. By including also the inputs, indicated by , the model reads


This relation is called the equation error model. On the right-hand side, the first term relating the past outputs at the present is referred to as the auto-regressive term; the second term reproduces the dynamics between the inputs and the outputs (exogenous part). The unknowns of the model are the set of coefficients , , the non-linear part – here represented by high-order polynomials – and the error . A classical way to model the error is to consider it as a moving-average of unknown white-noise. The most general strategy of identification based on Eq. 2 is the non-linear (N) autoregressive (AR) exogenous (X) with moving average (MA) algorithm, usually referred as N-ARMAX.

A quick overview of the algorithm is provided in A, while for a deeper discussion we refer to ljung1999system (); guidorzi2003multivariable (); semeraro2015techrep (). For our purposes, it is interesting to focus here on the choice of the inputs; in particular,

  1. if we are interested in revealing the core mechanisms of a dynamical system, we can neglect the input such that all the dynamics is described by the autoregressive part of the model and it is not contained in the inputs; in this case, the model is a discrete representation of an autonomous system.

  2. if a non-autonomous representation is sought, an optimal choice is represented by inputs modelled as stochastic white noise or random binary sequences. In this case, all the dynamics will be described by the identified input-output model.

  3. a third choice consists in introducing as input one of the available outputs , thus containing already some of the dynamics. In this case, due to the correlation between input and output, the resulting model will be an optimal transfer function among the elements; this approach is control-oriented.

In sec. 5, we use full non-linear ARMA models; the main goal is to identify the essential dynamics using autonomous models, in order to educe the intrinsic mechanisms at work. In sec. 6, we aim at identifying models that minimise the prediction error in terms of best mean-squared errors over a time-horizon, given a measurement as input ; in this case, full non-linear ARMAX models are computed.

The models discussed in the next sections were computed using the tools developed by Luis A. Aguirre and co-workers (see aguirre2009modeling ()). The linear models were cross-validated with in-house scripts semeraro2015techrep (). We use polynomial expansions, but alternative approaches are available. For instance, non-linear combinations of elementary functions can be used. In this case, the system identification problem changes: we do not seek only the coefficients of the model but also the optimal combination of functions fitting the original data; due to the necessity of identifying also the basis of functions, these techniques require full non-linear optimisation nelles2001nonlinear () or statistical-learning techniques (see SINDy (), for instance).

5 Non-linear models: dynamical system analysis from the observables

In this section we consider non-linear ARMA modelling. With respect to the full equation error model in Eq. 2, we neglect the exogenous terms. The basic model reads


The identification process consists in identifying the number of necessary terms and coefficients for an optimal representation of the time-series over a window of assimilation.

This approach leads to a discrete, autonomous model. The un-modelled part is accounted for in the error term . The final analysis is performed on a surrogate time-series obtained by running the model over long time-horizons, not driven from external known forcings. Thus, we do not expect a predictive model, but rather a phenomenological model. The aim is to mimic a part of the dynamics of the original system.

We analyse the results obtained by applying the identification approach to the time-series. In particular, we consider the TS-POD at , non-linearly filtered (see Sec. 3.1). As previously shown, the TS-POD series are less affected from measurement noise, in the range where the simultaneous data from TS-Series are available; moreover, the non-linear filter preserves the dynamics across all the frequencies. The chosen location corresponds roughly to the end of the potential core.

In Sec. 5.1, we provide a brief description of the parametric analysis performed during the computation of the non-linear ARMA models. Results are shown and discussed in Sec. 5.2 for three different models.

5.1 Computation of the chosen models and validation

The analysis involves several parameters, namely: i) the polynomial order; ii) the window of assimilation; iii) the maximum number of coefficients ; iv) the optimisation parameters. Due to the non-linear terms and the complexity of the system, we expect multiple minima for the optimisation process. In that sense, ranging a large span of parameters is necessary. In the following we briefly motivate and explain our choices.

Polynomial order – non-linearities

We consider as maximum non-linear order ; the choice is a compromise between numerical and computational issues rising from higher-order models and the physical description of the system. In particular, high-order non-linearities might lead to inconsistencies due to singular entries in the Hankel matrices (see A), used during the identification process.

Assimilation window

The model can depend from the chosen assimilation window. Indeed, considering the spectrogram in Fig. 6, we already observed that the dynamics is characterised by events at given frequencies corresponding to organised behaviour in time-domain. In principle, one should consider more windows of assimilation, stacked over the entire time-span; in practice, for computational and convergence reasons, we consider shorter time-windows of assimilation roughly corresponding to the average length of the high-amplitude events observed in the spectrogram. We tested windows of length ; in the present manuscript, we consider , that provided best results.


This coefficient dictates the second dimension of the regression matrix (see A), the first being the assimilation window. More importantly, it relates past entries with the actual value during the auto-regression process. We consider the range ; the maximum value in convective time corresponds to and it is based on the analysis of the auto-correlation. Larger values would lead to erroneous correlations between past and present dynamics.

Optimisation process parameters

For the optimisation process, we consider a maximum of iterations for each model and a maximum elements for the resulting models, whose terms are related to the error part.

As a result of the parametric analysis, we identified a total of models. As said, the models can be biased by the window of assimilation. At a first glance, one might think that the choice of relatively short time-windows of assimilation reduces the validity of the models; however, the main aim of the modelling approach is to explore the dynamics behind the organised behaviour observed in Fig. 6, trying to separate the inherent dynamics observable in small windows of observation from the effects of the noise-sensitivity. In that sense, it might not be possible to get endogenous models of the system, unless we were to introduce an exogenous forcing decorrelated from the output . For this reason, we believe that a short window of assimilation results in an identification process consistent with a part of the inherent dynamics of the dynamical system.

5.1.1 Validation process

The simulation of the models is not driven by known forcings; thus, these non-linear models will quickly deviate from the original time-series after a short transient roughly corresponding to the initial conditions. The validation process is performed in two steps: i) the minimisation of the error during the system identification; ii) the spectral analysis of the model, compared with the original time-series. In particular, we can not apply the Welch method due to the short time-window of assimilation; as an alternative, we apply a power spectral density estimation based on the Yule-Walker method stoica2005spectral (). This method estimates the spectral content by fitting an autoregressive linear prediction filter to the signal. By applying these two criteria and crossing the results, we were able to identify three different “families” of models. Of all the models identified, around show periodicity, quasi-periodicity, toroidal features or chaotic behaviour. The remaining models are badly conditioned, unstable or decaying to zero after a transient. The selected models are robust with respect to the initial conditions, i.e. they show the same asymptotic behaviour in the presence of uncertainties in the initial conditions. With respect to the position , the optimisation procedure leads to models with different coefficients , but similar dynamics from the phenomenological point of view. Thus, the surrogate models discussed here are still able to qualitatively describe the physics of the underlying system.

5.2 Autonomous models

In this section, we discuss three different models, labelled as follows

  1. Model T – in Fig. 11 the model is characterised by . It is a toroidal (T) example, that settles on the frequencies characterising the windows of learning.

  2. Model LC – in Fig. 12 the model is characterised by . After a rather long transient (up to ), it ends up on a stable limit-cycle (LC), see Fig. 13.

  3. Model S in Fig. 15 the model is characterised by . It is a chaotic model, as it can be seen from the first-return map. We believe that it reproduces some short (S) time dynamics related to the jittering.

The coefficient gives to the models distinctive features. The first two models are characterised by a long window of correlation. The third model is obtained considering the shortest, possible correlation length. In the following paragraphs, we comment on the significance of these phenomenological models.

(a) Comparison between time-series and model
(b) Autoregressive PSD
(c) Welch PSD
(d) Phase space comparison
(e) First-return map
Figure 11: Comparison N-ARMA model (in red for all the insets) vs. original time-series TS-POD at (in black). The comparison in time-domain is only qualitative. Thus, spectra analysis and phase space reconstruction are necessary for an assessment of the performance of the model. The embedding space for the phase space reconstruction is obtained with a delay . The first-return map is obtained using the maximum values in the model ; the black dots correspond to the beginning of the trajectory, while the superimposed blue dots the last points of the simulation.

5.2.1 Models and

In Fig. 11, the properties of the model are analysed. As already mentioned, none of these models is predictive: after a short time, the dynamics quickly deviate from the original one (Fig. 11). However, the original time-series are contaminated by the external disturbances, while the model behaviour is dictated by the intrinsic dynamic educed during the assimilation process. The quality of the model can be assessed by analysing the spectral properties in the window of assimilation: in Fig. 11, the power spectral density is shown based on autoregressive modelling; the current model reproduces the peak frequency ; a second peak appears at . However, we can also observe that model is not properly reproducing the low frequencies . The long-time behaviour is characterized by a toroidal behaviour in phase space, showed with a red solid line in Fig. 11. The Fourier analysis in Fig. 11 and the first-return map built with the relative maxima in Fig. 11 confirm the phase space representation. In particular, in Fig. 11, the black dots correspond to the beginning of the trajectory, while the superimposed blue dots the final part.

Model is characterised by a limited match in the autoregressive PSD; for this reason, a second model is presented in Fig. 12, labelled as model . In this case, we observe an excellent agreement at low Strouhal numbers, (see Fig. 12). We note a rather different behaviour when comparing the transient with the asymptotic state. During the transient, the spectrum is not characterized by distinct frequencies (see Fig. 12), although we observe the dominance of frequencies around and (due to the assimilation windows). The overall dynamics is more complex as shown in the phase space portrait in Fig. 12; the corresponding first-return map in Fig. 12 can be understood by considering the long-time behaviour of the system reported in Fig. 13. Indeed, after a transient of , the model settles on a limit cycle (Fig. 13); the dominating frequencies in Fig. 13 are the same as those observed in the Welch PSD analysed in Fig. 12. Having the long-time behaviour in mind, it is easier to understand the first-return map in Fig. 12. In particular, the blue dots are representative of the limit-cycle trajectory: we can observe that the model during the transient progressively settles on the limit cycle, that “attracts” the trajectory.

(a) Comparison between time-series and model
(b) Autoregressive PSD
(c) Welch PSD
(d) Phase space comparison
(e) First-return map
Figure 12: Comparison N-ARMA model (in red for all the insets) vs. original time-series TS-POD at (in black). The comparison in time-domain is only qualitative. Thus, spectra analysis and phase space reconstruction are necessary for an assessment of the performance of the model. The embedding space for the phase space reconstruction is obtained with a delay . The first-return map is obtained using the maximum values in the model ; the black dots correspond to the beginning of the trajectory, while the superimposed blue dots the last points of the simulation.
(a) Phase space
(b) First-return map
(c) Welch PSD
Figure 13: Asymptotic behaviour for N-ARMA model shown in Fig. 12: after the transient, the system settles on a stable limit cycle. In the phase space is shown (embedding delay ); in is shown the first-return map (blue dots of Fig. 12).
(a) at
(b) Model at
Figure 14: Spectrogram: comparison between the original time-series and the model .

The transient has a qualitative behaviour similar to the original model also if observed from the spectrogram. In Fig. 14, it is possible to compare the original time-series (inset ) with the model (), in the same time-span. It is evident that on average the length of the organised events is similar, while a significant difference is given by the frequency range; in particular, as already observed, this feature is an attribute of the extrinsic dynamics at work on the dynamical system. Thus, these two models suggest that if only the endogenous dynamics were considered, the dynamics due to non-linear interactions among wavepackets would be dominated by the organised behaviour that characterises the “short-living” events observed in the spectrogram. In particular, the dynamics would evolve after some transient on the surface of a -torus (T2) or a limit cycle. In models however, one limit cycle (at least) is stable and attracting. Therefore, the hypothetic underlying chaotic toroidal set is only transient.

5.2.2 Model S

In this section, we analyse a third model, depicted in Fig. 15 and based on a shorter window of correlation. Also in this case, the peak in frequency attains at , as shown in the auto-regressive PSD analysis in Fig. 15; this same dominance is observed in the Welch analysis in Fig. 15. The resulting behaviour in phase space and in the first-return map is rather distinctive and different compared to the previous models. In particular, this model appears to slowly evolve in time; in this case, we can conjecture that the model captures some inherent mechanisms related to combination of wavepackets-pairs being driven by the stochastic turbulence in the real case.

5.2.3 Discussion

Within the limits of the methodology, an analysis of the models provides some hints regarding the possible dynamics underlying the system. The main idea is that the jittering is a mixture of stochastic and deterministic components. Model T and Model LC seems to indicate that frequencies can lock together, leading to non-linear wave-wave interactions. Eventually, without external forcing, they can set on limit cycles or in trajectories embedded on a toroidal surface in phase space.

(a) Comparison between time-series and model
(b) Autoregressive PSD
(c) Welch PSD
(d) Phase space comparison
(e) First-return map
Figure 15: Comparison N-ARMA model (in red for all the insets) vs. original time-series TS-POD at (in black). The comparison in time-domain is only qualitative. Thus, spectra analysis and phase space reconstruction are necessary for an assessment of the performance of the model. The embedding space for the phase space reconstruction is obtained with a delay . The first-return map is obtained using the maximum values in the model ; the black dots correspond to the beginning of the trajectory, while the superimposed blue dots the last points of the simulation.

In this sense, the jitter signature is now more deterministic, but, in the real case, it is ’randomised’ due to continuous forcing by stochastic background turbulence. In fact, because the flow is a convective amplifier, the wave-wave interactions can never fully take over and the asymptotic behaviour can never happen because information is continuously convected out of the system. The dynamical systems reproduced by Model T and Model LC evolve according to initial conditions and a dynamic law that cannot evacuate the information convectively.

One might imagine the phase space as a dense ensemble of trajectories evolving on limit cycles and toroidal surfaces, associated with wave-wave interactions. Due to the sensitivity of the system to external noise, the system will continuously and erratically deviate from these attracting solutions without setting on any one of them as observed in the spatiotemporal plots and the spectrogram in Sec. 3. The high sensitivity to external disturbances is inherently connected to the convective nature of the flow. This picture can be linked to recent ideas regarding wavepacket theory: the organisation observed is understood as being due to the resolvent operator forced by background turbulence. In this sense, the resolvent operator is an organiser of the flow, while the turbulence that forces wavepackets via this operator is the disorganiser.

Also, note that the resolvent analysis has proved to be particularly efficient for , where the linear wavepackets are in remarkable agreement with the coherent structures observable in the jet flow. On the other hand, at low – and particularly at low Mach number – non-linear wave-wave interactions may be more relevant as the linear analogy fails: this makes our conjecture particularly appealing for explaining this discrepancy.

6 Input-output modelling for control: is linear modelling enough?

When the forcing term is represented by an output of the system, used as an input of the model, system identification identifies a transfer function between inputs and outputs. In particular, we use the output at and inputs placed in four positions, . The models include only one input (single-input-single-output modelling, SISO); the performance are shown in Fig. 16-, where we compare the model prediction with the original time-series. The window of assimilation is . The validation is performed by integrating the model in the window .

The application of the full non-linear ARMAX algorithm results in a polynomial expression: the final models are fully linear and based on a limited number of coefficients and , . The good results can be explained by observing the spectral analysis in Fig. 5. As already noted, the convective nature of the flow implies a progressive drift of the main frequency of amplification and a slight variation of the maximum amplitude in the spectrum as we move further downstream; in this sense, we need only a limited number of degrees of freedom for modulating amplitudes and frequencies, and correctly representing the time-delays of the system. This observation justifies the small dimensions of the models, that here act as filters.

The standard deviation between the resulting model and the original time-series is generally low. In Fig. 17, the results are summarised for a number of outputs along the streamwise direction. In particular, six different positions along are introduced as input, in the range , while as outputs are considered the locations downstream of the inputs between and . Except the transfer functions with inputs at and output at , characterised by a standard deviation , in all the other cases the standard deviation is below . As already mentioned above, all the shown models are linear: this observation has interesting implications when considering the design of linear controllers.

(a) Input at
(b) Input at
(c) Input at
(d) Input at
Figure 16: Reduced-order model vs. original time-series. Solid, black lines indicate the original signal at ; the red dots indicate the model prediction. Four inputs are considered, at .

These results confirm and extend the findings of sasaki2015real (), where linear ARMAX modelling, parabolised stability equation (PSE) and empirical transfer function were compared for the estimation of downstream propagating wavepackets, showing equivalent results. With respect to the work of sasaki2015real (), in our analysis the linearity of the model is a result of the optimisation process and not a working hypothesis: despite the fact that candidate models introduced in the identification procedure were characterised by all the terms of Eq. 2 up to order , non-linearities appear to be un-necessary for an optimal representation of the transfer functions.

6.1 Robustness and control

The limited amount of dynamical information retained in these models lead to robustness issues. The set of coefficients is obtained as least-squares solution from the assimilation datasets at given amplitude. The transfer function only retains information regarding the time-delays of the system, the modulation of the amplitude and the frequency bandwidth. Any sources of uncertainty can lead to robustness issues. As an example, in Fig. 17, we illustrate the lack of robustness with respect to time-delays, due to the convective time of propagation of the perturbations. The model assimilated using the input at is applied with different inputs, taken in the range . The validity of the models quickly decays. Values of the standard deviation indicate phase-opposition of the resulting model with respect to the original time-series.

The robustness issues can be attacked using different strategies. Control design can be conceived by accounting for the robustness issues arising from uncertainties of the system within the robust framework or the adaptive control theory (see sipp2016linear (); fabbiane2014adaptive ()). In particular, with respect to the convective nature of the flow, feed-forward schemes combined with self-tuning, adaptive controllers can be applied. These controllers are characterised by two time-scales: a fast time-scale, for the real-time estimation of the input , and a slow time-scale, based on feedback measurements of the flow necessary for correcting off-design conditions. Examples of these strategies can be found for weakly non-linear transitional flows (see semeraro2013transition ()) or non-linear, turbulent cases (see rathnasingham2003active ()).

(a) Models’ performance
(b) Robustness
Figure 17: Left: the normalized standard deviation between the models and the time-series is analysed. The -axis indicates the location of the output, while in the legend the six lines correspond to the inputs. In total, models are analysed, most of which are characterised by a standard deviation . Right: the model is computed between the input at and the output at ; the same model, is tested by using different inputs between and ; the test shows the lack of robustness with respect to the time-delays of the system.

7 Conclusions

Coherent large-scale structures in turbulent jets – usually referred as wavepackets – have been subject of numerous investigations for the role that such structures might play in sound radiation. Remarkably, it is well established that certain aspects of the wavepacket dynamics can be modelled by means of linear (jordan2013wave ()). In this work, we considered an unforced, isothermal, subsonic jet issuing from a round nozzle with a fully turbulent boundary layers at and . The temporal behaviour of the axisymmetric wavepackets is analysed; starting from the studies by breakey2013near (), we investigate the experimental, spatio-temporal measurements from a dynamical point of view. To the best of our knowledge, this is the first time such an analysis has been attempted. The key point is the remarkable spatio-temporal coherence of the measurements: despite the high Reynolds number, the low-dimensional projection of the dynamics reveals a well-organised behaviour. This aspect raises questions regarding the temporal behaviour of the wavepackets, the nature of the non-linearities involved and the possibility of designing control strategies starting from this knowledge.

In order to answer these questions, a variety of strategies ranging from statistical tools to system identification have been adopted. The estimation of the correlation dimension and the embedding dimension of the dynamical system confirmed that the organised temporal behaviour is low-dimensional and deterministic. The minimal embedding ranges between and , while the correlation dimension is rather low at locations corresponding to the end of the potential core, . These dimensions are nonetheless still too large for useful characterisation of attractor geometry.

For this reason, alternative strategies are necessary for exploring the dynamics of the system. A system identification approach is pursued. The time-series taken at a streamwise location of the jet is assumed as output; we replace the original dynamics with surrogate, autonomous models using data-assimilation based on this output. These models are limited, as they reproduce only a part of the dynamics and are influenced by the choice of the parameters.

By comparing the results from three different classes of models, we suggest that the jittering dynamics might be described as a combination of extrinsic and intrinsic mechanisms at work. We observe that our surrogate models reproduce the ideal case where the background turbulence is not active; in this limit, the dynamics of the system would be dominated by wave-wave interactions leading to stable limit cycles and/or solutions lying on toroidal attractors. Evidence on this behaviour is provided by the spectral analysis of the transient dynamics of these solutions that qualitatively reproduces the local (in time) dynamics of the original system. However, in the real case, the extrinsic dynamics drives the system: wave-wave interactions can never fully take over and the asymptotic behaviour of the ideal models can never happen. This interpretation justifies the self-organised features appearing in the time-frequency analysis of the time-series and it is consistent with the nature of the jets as noise-amplifiers: due to the sensitivity of the system to external noise, the trajectory in phase space will continuously escape from these attracting solutions.

Given these interpretations, control strategies can be conceived. Due to the convective nature of the system and its implications regarding the non-linear interactions, system identification is used for computing non-autonomous models, where the inputs are represented by local measurements taken upstream of the outputs. These models are filters that i) account for the time-delays and ii) modulate the amplitudes and the frequencies already existing in the inputs. In particular, we found that non-linearities are un-necessary for an optimal representation of these models. In other words, if the inputs contain a large part of the dynamics, a linear filter is enough to reproduce the temporal behaviour at downstream locations due to the convective nature of the system. The main limit is the robustness of these linear filters; however, in our opinion, this might be more a problem of controllability of the flow, rather than a lack observability and predictability. Thus, we believe that the robustness limitations may be tackled in the control design, by applying strategies such as adaptive control or loop-shaping robust control, and verifying the effects of the actuation on the real flow.

Future work should be devoted to further improving the non-linear models for the dynamical analysis. While accepting the impossibility of capturing all the features of the system in low-dimensional models, a possible path to improve our understanding of the jittering would be the identification of surrogate models where a stochastic excitation is also included; in principle, this strategy would help in separating the instrinsic mechanisms from the extrinsic ones. In this sense, we believe that machine learning and statistical learning techniques may be robust alternatives to classical non-linear system identification. Indeed, these techniques are quickly spreading in the community as a powerful tool for system identification (see SINDy ()).

The authors wish to acknowledge Luis A. Aguirre for sharing his system identification package and Christophe Letellier for fruitful discussions. This work is supported by the Agence Nationale de la Recherche (ANR) under the ”Cool Jazz” project, grant number ANR-12-BS09-0024.

Appendix A System identification using polynomials expansion

We introduced in Sec. 4 the basic ideas behind the application of system identification for the qualitative analysis of dynamical system and the definition of control-oriented models. The aim of the section is to provide a quick overview of the algorithms.

We rely on an input-output description of time-invariant systems. The aim of our models is reproducing the behaviour of the output by defining a non-linear model of the dynamics of the measurement or a transfer function between the input and the output . For sake of conciseness, we consider a single-input-single-output (SISO) case; nonetheless, the extension to the multi-inputs-multi-outputs cases is rather straightforward from the theoretical point of view.

By assuming a SISO system, we can generalise the equation error models in terms of Volterra series


where are the order Volterra kernel. The error can be modelled as a moving-average of (unknown) white-noise


where are the coefficients for the error modelling. The kernel’s coefficients and the error are the unknowns of our problem. The Volterra-series are an extension of the linear convolution integral and represent non-linear systems as a series of multi-summations (or integrals) of the Volterra kernels and the inputs billings2013nonlinear (). We can note that the complete model consists of a linear combination of three groups of terms

  1. the auto-regressive terms related to the output ;

  2. the exogenous terms due to the input ;

  3. the moving-average description of the error .

Additionally, the non-linear combinations are introduced in the full model.

a.1 Solution of the identification problem: a basic ARMAX algorithm

The first step of the identification procedure consists in a the preliminary computation of the coefficients . To this end, we assume that the error is null, such that the least-squares regression can be performed. Note that the error is modelled as a moving average of white noise , leading to a coloured noise with expected value ; the estimate of the coefficients will be affected by a bias, due to the expected value of the error . In particular, the error is computed as


where is a rectangular Hankel matrix including inputs and outputs . The Hankel matrix has the dimensions of the number of the coefficients and the window of assimilation. We can introduce a matrix of the same dimensions as the Hankel matrix such that


Thus, with a proper choice , the “de-biased” estimate of can be computed as


if is invertible. This generalisation of the least-squares problem is known as Instrumental Variable (IV) method and allows a first estimation of the coefficients in the expansion, while neglecting the error.

The IV method is iterative. For the first iteration, the matrix of instruments is written as ; a rather effective choice consists in choosing , i.e. the delayed input. These values are relaxed till convergence.

The IV method allows to compute the coefficients and , reducing the effects of the bias due to the non-whiteness of the error . The unknowns of the moving-average part, namely the coefficients and the white-noise series can be extracted using appropriate algorithms (see guidorzi2003multivariable () or semeraro2015techrep ()).

Thus, we can divide the basic algorithm into two steps:

  1. Evaluate the coefficients and using the iterative IV method.

  2. Apply the moving-average algorithm for defining the coefficients and the white-noise .

This prototypical algorithm allows a first estimate of the kernel coefficients and the error of the moving average description. We can improve these estimates using the extended least-squares (ELS) approach; in particular, we can consider the first estimation of the error as an input of the system: in this way, the estimates will be de-biased by explicitly accounting the error. Moreover, an iteration can be cast: the difference between the prediction and the original output, allows to refine again the kernels estimate and further improve the prediction.

a.2 Other techniques

The main difference between the linear modelling and the full non-linear Volterra series is mainly related to the dimensions of the problem and the risk of ill-conditioness. Already cubic non-linearities may introduce nearly null entries, leading to rank-deficiency and singularities. To this end, the regression step can be replaced by well-conditioned algorithms, such as the orthogonal least-squares (OLS); this iterative scheme introduces in the least-squares process a Gram-Schmidt orthogonalisation, which ensures that each new column added in the matrices is orthogonal to all previous columns. Moreover, this process allows us to estimate the relative relevance of the kernel’s coefficient in residual terms; thus, elements that are not relevant are discarded during the process. The OLS approach is described by billings2013nonlinear () or aguirre2009modeling (), and allows the identification of minimal non-linear models by keeping only the relevant terms of the model. This is the strategy implemented in the routines adopted in our manuscript.

Other methods for the system identification are based on the minimisation of the prediction at a given time-horizon, based on gradient-based optimisation process; an example is given by the predictive error method (PEM), implemented in the system identification toolbox by ljung2007system () (see also ljung1999system ()).



  • (1) P. Jordan, T. Colonius, Wave packets and turbulent jet noise, Annual Review of Fluid Mechanics 45 (2013) 173–195.
  • (2) D. E. S. Breakey, P. Jordan, A. V. G. Cavalieri, O. Léon, M. Zhang, G. Lehnasch, T. Colonius, D. Rodriguez, Near-field wavepackets and the far-field sound of a subsonic jet, in: 19th AIAA/CEAS Aeroacoustics Conference, AIAA Paper, Vol. 2083, 2013, pp. –.
  • (3) E. Mollö-Christensen, R. Narasimha, Sound emission from jets at high subsonic velocities, Journal of Fluid Mechanics 8 (01) (1960) 49–60.
  • (4) E. Mollö-Christensen, M. A. Kolpin, J. R. Martucelli, Experiments on jet flows and jet noise far-field spectra and directivity patterns, Journal of fluid mechanics 18 (02) (1964) 285–301.
  • (5) E. Mollo-Christensen, Jet noise and shear flow instability seen from an experimenter’s viewpoint, Journal of Applied Mechanics 34 (1967) 1–7.
  • (6) S. C. Crow, F. H. Champagne, Orderly structure in jet turbulence, Journal of Fluid Mechanics 48 (03) (1971) (547–591).
  • (7) J. C. Lau, M. J. Fisher, H. V. Fuchs, The intrinsic structure of turbulent jets, Journal of Sound and Vibration 22 (04) (1972) 379–384.
  • (8) A. Michalke, H. V. Fuchs, On turbulence and noise of an axisymmetric shear flow, Journal of Fluid Mechanics 70 (1975) 179–205.
  • (9) R. R. Armstrong, A. Michalke, H. Fuchs, Coherent structures in jet turbulence and noise, AIAA Journal 15 (7) (1977) 1011–1017.
  • (10) C. J. Moore, The role of shear-layer instability waves in jet exhaust moise, Journal of fluid mechanics 80 (1977) (321–367).
  • (11) A. Michalke, Instability of a compressible circular free jet with consideration of the influence of the jet boundary layer thickness, Z. fur Flugwissenschaften,(West-Germany) 19 (8) (1971) (319–328).
  • (12) D. Bechert, E. Pfizenmaier, On the amplification of broad band jet noise by a pure tone excitation, Journal of Sound and Vibration 43 (3) (1975) 581–587.
  • (13) J. Cohen, I. Wygnanski, The evolution of instabilities in the axisymmetric jet. part 1. the linear growth of disturbances near the nozzle, Journal of Fluid Mechanics 176 (1987) 191–219.
  • (14) D. G. Crighton, M. Gaster, Stability of slowly diverging jet flow, Journal of Fluid Mechanics 77 (02) (1976) (397–413).
  • (15) P. Plaschko, Helical instabilities of slowly divergent jets, Journal of Fluid Mechanics 92 (02) (1979) 209–215.
  • (16) A. Michalke, A wave model for sound generation in circular jets, Zentralstelle für Luftfahrtdokumentation und-information, 1970.
  • (17) A. Michalke, An expansion scheme for the noise from circular jets, Tech. rep., DTIC Document (1971).
  • (18) K. A. Bishop, J. E. Ffowcs-Williams, W. Smith, On the noise sources of the unsuppressed high-speed jet, Journal of Fluid Mechanics 50 (01) (1971) 21–31.
  • (19) S. C. Crow, Acoustic gain of a turbulent jet, Phys. Soc. Meeting, Univ. Colorado, Boulder, paper IE 6.
  • (20) D. G. Crighton, Basic principles of aerodynamic noise generation, Prog. Aerospace Sci. 16.
  • (21) D. G. Crighton, P. Huerre, Shear-layer pressure fluctuations and superdirective acoustic sources, Journal of Fluid Mechanics 220 (1) (1990) 355–368.
  • (22) T. Suzuki, T. Colonius, Instability waves in a subsonic round jet detected using a near-field phased microphone array, Journal of Fluid Mechanics 565 (2006) 197–226.
  • (23) A. V. G. Cavalieri, D. Rodríguez, P. Jordan, T. Colonius, Y. Gervais, Wavepackets in the velocity field of turbulent jets, Journal of Fluid Mechanics 730 (2013) 559–592.
  • (24) P. Jordan, T. Colonius, G. A. Bres, M. Zhang, A. Towne, S. K. Lele, Modeling intermittent wavepackets and their radiated sound in a turbulent jet, in: Proceedings of the Summer Program, 2014, p. 241.
  • (25) J. W. Nichols, S. K. Lele, Global modes and transient response of a cold supersonic jet, Journal of Fluid Mechanics 669 (2011) 225–241.
  • (26) X. Garnaud, L. Lesshafft, P. Schmid, P. Huerre, The preferred mode of incompressible jets: linear frequency response analysis, Journal of Fluid Mechanics 716 (2013) 189–202.
  • (27) O. Schmidt, A. Towne, T. Colonius, P. Jordan, V. Jaunet, A. V. Cavalieri, G. A. Brès, Super-and multi-directive acoustic radiation by linear global modes of a turbulent jet, in: 22nd AIAA/CEAS Aeroacoustics Conference, 2016, p. 2808.
  • (28) G. A. Brès, J. Jaunet, M. Le Rallic, P. Jordan, T. Colonius, S. K. Lele, Large eddy simulation for jet noise: the importance of getting the boundary layer right, in: 21st AIAA/CEAS Aeroacoustics Conference, 2015, p. 2535.
  • (29) M. Zhang, P. Jordan, G. Lehnasch, A. V. G. Cavalieri, A. Agarwal, Just enough jitter for jet noise, in: 20th AIAA/CEAS Aeroacoustics Conference, AIAA Paper, Vol. 3061, 2014, pp. –.
  • (30) Y. B. Baqui, A. Agarwal, A. V. Cavalieri, S. Sinayoko, A coherence-matched linear source mechanism for subsonic jet noise, Journal of Fluid Mechanics 776 (2015) 235–267.
  • (31) A. V. G. Cavalieri, P. Jordan, A. Agarwal, Y. Gervais, Jittering wave-packet models for subsonic jet noise, Journal of Sound and Vibration 330 (18) (2011) 4474–4492.
  • (32) A. V. G. Cavalieri, G. Daviller, P. Comte, P. Jordan, G. Tadmor, Y. Gervais, Using large eddy simulation to explore sound-source mechanisms in jets, Journal of Sound and Vibration 330 (17) (2011) 4098–4113.
  • (33) F. Kerhervé, P. Jordan, A. V. G. Cavalieri, J. Delville, C. Bogey, D. Juvé, Educing the source mechanism associated with downstream radiation in subsonic jets, Journal of Fluid Mechanics 710 (2012) 606–640.
  • (34) A. Towne, T. Colonius, P. Jordan, A. V. G. Cavalieri, G. A. Brès, Stochastic and nonlinear forcing of wavepackets in a mach 0.9 jet, in: 21st AIAA/CEAS Aeroacoustics Conference, 2015, p. 2217.
  • (35) G. Tissot, M. Zhang, F. Lajus Jr., A. V. G. Cavalieri, P. Jordan, T. Colonius, Sensitivity of wavepackets in jets to non-linear effects: the role of the critical layer, Journal of Fluid Mechanics In review.
  • (36) N. D. Sandham, C. L. Morfey, Z. W. Hu, Nonlinear mechanisms of sound generation in a perturbed parallel jet flow, Journal of Fluid Mechanics 565 (2006) 1–23.
  • (37) V. Suponitsky, N. D. Sandham, C. L. Morfey, Linear and nonlinear mechanisms of sound radiation by instability waves in subsonic jets, Journal of Fluid Mechanics 658 (2010) 509–538.
  • (38) K. Sasaki, S. Piantanida, A. Cavalieri, P. Jordan, Real-time modelling of wavepackets in turbulent jets, in: 21nd AIAA/CEAS Aeroacoustics Conference, AIAA Paper, 2015.
  • (39) K. Sasaki, G. Tissot, A. Cavalieri, S. Silvestre, P. Jordan, D. Biau, Closed-loop control of wavepackets in a free shear flow, in: 22nd AIAA/CEAS Aeroacoustics Conference, AIAA Paper, Vol. In review, 2016.
  • (40) N. Gautier, J.-L. Aider, T. Duriez, B. Noack, M. Segond, M. Abel, Closed-loop separation control using machine learning, Journal of Fluid Mechanics 770 (2015) 442–457.
  • (41) A. V. G. Cavalieri, P. Jordan, T. Colonius, Y. Gervais, Axisymmetric superdirectivity in subsonic jets, Journal of fluid Mechanics 704 (2012) 388–420.
  • (42) S. Narayanan, G. H. Gunaratne, F. Hussain, A dynamical systems approach to the control of chaotic dynamics in a spatiotemporal jet flow, Chaos: An Interdisciplinary Journal of Nonlinear Science 23 (3) (2013) 033133.
  • (43) F. Takens, Detecting strange attractors in turbulence, in: Dynamical systems and turbulence, Warwick 1980, Springer, 1981, pp. 366–381.
  • (44) P. Grassberger, I. Procaccia, Characterization of strange attractors, Physical review letters 50 (5) (1983) 346.
  • (45) H. Kantz, T. Schreiber, Nonlinear time series analysis, Vol. 7, Cambridge university press, 2004.
  • (46) L. Cao, Practical method for determining the minimum embedding dimension of a scalar time series, Physica D: Nonlinear Phenomena 110 (1) (1997) 43–50.
  • (47) R. Hegger, H. Kantz, T. Schreiber, Practical implementation of nonlinear time series methods: The tisean package, Chaos: An Interdisciplinary Journal of Nonlinear Science 9 (2) (1999) 413–435.
  • (48) J.-P. Eckmann, D. Ruelle, Fundamental limitations for estimating dimensions and lyapunov exponents in dynamical systems, Physica D: Nonlinear Phenomena 56 (2-3) (1992) 185–187.
  • (49) P. Huerre, P. A. Monkewitz, Local and global instabilities in spatially developing flows, Annual review of fluid mechanics 22 (1) (1990) 473–537.
  • (50) L. A. Aguirre, C. Letellier, Modeling nonlinear dynamics and chaos: a review, Mathematical Problems in Engineering 2009.
  • (51) L. Ljung, System identification, Wiley Online Library, 1999.
  • (52) R. Guidorzi, Multivariable system identification: from observations to models, Bononia University Press Bologna, 2003.
  • (53) O. Semeraro, L. Mathelin, An open-source toolbox for linear system-identification., Tech. rep., LIMSI-CNRS, Orsay, France (2016).
  • (54) O. Nelles, Nonlinear system identification: from classical approaches to neural networks and fuzzy models, Springer Science & Business Media, 2001.
  • (55) S. L. Brunton, J. L. P. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113 (15) (2016) 3932–3937.
  • (56) P. Stoica, R. L. Moses, Spectral analysis of signals, Vol. 452, Pearson Prentice Hall Upper Saddle River, NJ, 2005.
  • (57) D. Sipp, P. J. Schmid, Linear closed-loop control of fluid instabilities and noise-induced perturbations: A review of approaches and tools, Applied Mechanics Reviews 68 (2) (2016) 020801.
  • (58) N. Fabbiane, O. Semeraro, S. Bagheri, D. S. Henningson, Adaptive and model-based control theory applied to convectively unstable flows, Applied Mechanics Reviews 66 (6) (2014) 060801.
  • (59) O. Semeraro, S. Bagheri, L. Brandt, D. S. Henningson, Transition delay in a boundary layer flow using active control, J. Fluid Mech. 731 (9) (2013) 288–311.
  • (60) R. Rathnasingham, K. S. Breuer, Active control of turbulent boundary layers, Journal of Fluid Mechanics 495 (2003) 209–233.
  • (61) S. A. Billings, Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains, John Wiley & Sons, 2013.
  • (62) L. Ljung, System identification toolbox for use with MATLAB.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description