Mean- and unsteady-flow reconstruction with one or two time-resolved measurements
MMean- and unsteady-flow reconstruction with one or twotime-resolved measurements
Lucas Franceschini a , Denis Sipp a , Olivier Marquet a a Department of Aerodynamics, Aeroacoustics and Aeroelasticity, ONERA, Meudon, France
February 9, 2021
Abstract
In this article, we propose a methodology to reconstruct, in a single step, the mean- and unsteadyproperties of a flow from very few time-resolved measurements. The procedure is based on the apriori alignement of Fourier- and Resolvent-modes over energetic frequencies, which is a commonfeature in shear-dominated transitional flows. Hence, the Reynolds-stresses, which determine themean-flow, may be approximated from a series of Resolvent modes, which discretize the fluctuationfield in the frequency space and whose amplitudes can be tuned thanks to few measurements. Inpractice, we solve a nonlinear optimization problem (with only few parameters) based on a modelcoupling strongly the equations governing the mean-flow and the Resolvent modes. The input datafor the assimilation procedure may be very sparse, typically one or two pointwise measurements.This technique is applied to two distinct physical configurations, one ”oscillator” flow with periodicfluctuations (squared-section cylinder) and one ”noise amplifier” flow with a broadband frequencyspectrum (backward-facing step).
Keywords:
Data Assimilation, Stability Analysis, Mean-flow Analysis
The reconstruction of a time-averaged (mean) flow and its fluctuations is a topical issue within data-assimilation influid mechanics. Mean-flow reconstruction alone can be performed through a variational minimization algorithm,where tuning parameters of a mean-flow model are adjusted to match available measurement data from experi-ments or higher-fidelity numerical simulations. Foures, Dovetta, Sipp, and Schmid (2014) applied this techniqueon the incompressible laminar flow over a circular cylinder. They used a spatially-distributed volume forcingterm (modeling the divergence of the Reynolds-stress tensor) as control parameter to reconstruct the laminarmean-flow using mean-velocity measurements. The same technique was employed by Symon, Dovetta, McK-eon, Sipp, and Schmid (2017) on an experimental incompressible flow around an idealized airfoil at much higherReynolds number. One drawback of this approach is the need for a large number of well-placed measurementsfor the reconstruction to produce accurate results. This issue was partially addressed by Foures et al. (2014) whostudied the ability of the algorithm to interpolate/extrapolate the field in regions of the physical domain far fromthe measurements. In the case of high-Reynolds number turbulent flows, it becomes mandatory to consider aturbulence model. Such a strategy was initially introduced for the field-inversion step to learn new turbulencemodels (Duraisamy, Iaccarino, & Xiao, 2019): for example, Ferrero, Iollo, and Larocca (2020); He, Liu, Gan, andLesshafft (2019); Parish and Duraisamy (2017); Singh and Duraisamy (2016) considered a spatially distributedcorrection term acting on the turbulence production in the Spalart-Allmaras model. Franceschini, Sipp, and Mar-quet (2020) actually showed that it may be better to consider a force correction term in the momentum equationswhen dense measurements are available and a scalar source term in the turbulence equation in the case of sparsemeasurements.In the case of dense measurements and a force correction term, Foures et al. (2014); Franceschini etal. (2020); Symon et al. (2017) showed that the procedure may lead to an accurate reconstruction of the Reynoldsstress tensor (from the force term).Those techniques remain purely steady approaches, where mean-flow information is used to reconstruct themean-flow field. If we target the reconstruction of the fluctuating field in the case of few pointwise measurements,we need to resort to unsteady data-assimilation techniques. Classically, those methods belong to two (possiblyoverlapping) categories: 4DVar (Le Dimet & Talagrand, 1986) and Ensemble Kalman filters (Evensen, 2009).Those methods have been extensively applied in the context of meteorology (Liu, Xiao, & Wang, 2008; Lorenc,1986) and more recently in fluid mechanics (A., D., & E., 2013; D’Adamo, Papadakis, Memin, & Artana, 2007;Meldi, 2018; Mons, Chassaing, Gomez, & Sagaut, 1986). Such techniques are computationally intensive sincethey are based on an unsteady model that is marched forward in time, while the optimization step is performed a r X i v : . [ phy s i c s . f l u - dyn ] F e b ither by a time-adjoint code (for 4DVar) and/or a statistical treatment of an ensemble of points in the phasespace (for Ensemble Kalman filter).An alternative approach for the reconstruction of the fluctuation field relies on the analysis of the linearNavier-Stokes operator around the mean-flow. The origin of this approach lies in the works of Barkley (2006);Pier (2002), who observed that, in the case of supercritical cylinder flow, the stability analysis of the linearizedNavier-Stokes operator on the time-averaged (mean) flow (instead of the steady solution or the base-flow) givesrise to a mode whose eigenvalue is purely imaginary and its frequency close to the one of the nonlinear limit cycle.The theoretical explanation of this fact was further explored by Sipp and Lebedev (2007) with the aid of weakly-nonlinear analysis close to the bifurcation and then by Mezi´c (2013); Turton, Tuckerman, and Barkley (2015),who showed that the Fourier mode of a harmonic (in contrast to time-periodic) flow corresponds to a marginaleigenmode of the Jacobian operator. Later, this led Mantic-Lugo, Arratia, and Gallaire (2015) to develop aself-consistent model describing the nonlinear saturation of such oscillator flow, where the steady solution evolvestoward the mean-flow through the action of the Reynolds-stresses, which are approximated from the marginaleigenmode associated to the mean-flow. The extension of the mean-flow analysis to flows characterized by abroadband spectral frequency content has been conducted by Hwang and Cossu (2010); McKeon and Sharma(2010) in streamwise homogeneous flows and by Beneddine, Sipp, Arnault, Dandois, and Lesshafft (2016) inmore general configurations. They showed that, under some fairly general assumptions (in shear dominatedflows, strong spatial convective instabilities such as the Kelvin-Helmholtz instability favour rank-1 Resolventoperators), the Fourier mode of the nonlinear fluctuation field at a given frequency is proportional to the leadingmode of the Singular Value Decomposition of the Resolvent operator around the mean-flow, evaluated at thatfrequency. In stochastic flows, Towne, Schmidt, and Colonius (2018) actually showed that the Resolvent modesare equal to the Spectral Proper Orthogonal Decomposition (SPOD) modes in the case of a white noise spatiallyuncorrelated forcing. Yet, this is generally not the case in turbulent flows and the alignement between Resolventand SPOD modes may be less perfect. The use of an eddy-viscosity field modeling the small-scale turbulencein the Resolvent analysis may strongly increase the alignement of these modes (Morra, Semeraro, Henningson,& Cossu, 2019): Pickering, Rigas, Schmidt, Sipp, and Colonius (2020) showed that it is possible to use andadjoint-based optimization strategy to determine a single eddy-viscosity field that simultaneously improves thealignement of several Resolvent modes at a given frequency in highly turbulent compressible jet flows.Coming back to the estimation problem, if the Fourier mode and Resolvent modes are assumed to be alignedat a given frequency, we know the spatial structure of the fluctuation field at that frequency, up to a (complex)multiplicative constant which represents the energy and phase of that mode. Beneddine et al. (2016); Beneddine,Yegavian, Sipp, and Leclaire (2017); G´omez, Blackburn, Rudman, Sharma, and McKeon (2016); Towne, Lozano-Dur´an, and Yang (2020) showed how to estimate this constant with few point-wise time-resolved measurements.From a practical viewpoint, the main drawback of this approach is of course the need for the a priori knowledgeof the mean-flow. He et al. (2019); Symon, Sipp, and McKeon (2019); Symon, Sipp, Schmid, and McKeon (2020)therefore proposed to combine in serial a variational mean-flow data-assimilation to determine an approximatemean-flow and a Resolvent analysis for the estimation of the unsteady fluctuation field.In the present work, we propose to combine in a single step both the estimation of the mean-flow and theestimation (by Resolvent analysis) of the fluctuation field. By doing so, we obtain a nonlinear model whoseunknown parameters are the amplitudes of the Resolvent modes at each considered frequency. To tune theseparameters, we only need few sensors (of the same order as the number of unknowns). For this reason, one couldconsider this model as a reduced order model for the unsteady flow, which is then used for data-assimilation. Onlyvery sparse information is required here for the reconstruction of the mean- and fluctuating quantities. This is astrong improvement with respect to the uncoupled serial approach, since an accurate recovery of the mean-flowrequires denser measurements (at least lines of measurements). This procedure will be applied to simple 2Dmodel problems of transitional laminar flows, leaving therefore aside the need of an eddy-viscosity and the issueof multiple strong SPOD modes: we assume as in Beneddine et al. (2017) that the alignment between dominantResolvent modes and Fourier modes (if only one SPOD mode is energetic at a given frequency, we will call it aFourier mode) is high for all energetic frequencies). The first configuration will be the flow around a square cylinder( Re = 100), falling in the category of oscillator flows, for which the frequency spectrum exhibits a dominant peakand smaller harmonics. This means that only one Resolvent mode may be sufficient at the fundamental frequency(and thus only one parameter) for the reconstruction procedure. We will show that the velocity at some locationof the flow then suffices for the estimation of the mean-flow and the fluctuation field. The second configurationis the two-dimensional backward-facing step ( Re = 500), a typical example of noise amplifier Herv´e, Sipp, andSchmid (2012), where unsteadiness is triggered by upstream located white-noise in time forcing. Here again, thefrequency content of the flow at some points can be used to tune a series of amplitudes of Resolvent modes, whichdiscretize the broadband fluctuation field in the frequency space.The article proceeds as follows. First ( § Theory
This section is devoted to the derivation of the model that will be used in the data-assimilation procedure. Weconsider an incompressible flow, whose velocity and pressure fields are governed by the Navier-Stokes equationswritten in compact form: ∂ t u + N ( u ) = , where N ( u ) = ( u · ∇ ) u + ∇ p − Re − ∆ u (1)is the nonlinear Navier-Stokes operator, with the pressure p such that the velocity field is incompressible ∇· u = 0.The Reynolds number Re = U ∞ L/ν is based on an arbitrary velocity U ∞ and a characteristic length L , which areused to make all the variables non-dimensional. Since we focus on statistically steady regimes, we may decomposethe instantaneous state into a mean and a fluctuating component, as u ( x , t ) = u ( x ) + u (cid:48) ( x , t ) , p ( x , t ) = p ( x ) + p (cid:48) ( x , t ) , (2)where the time-average operator is defined as( · ) = lim T → + ∞ T (cid:90) T ( · )( t ) dt. (3)Introducing the decomposition (2) into the Navier-Stokes equations (1) and taking the time-average yields themean flow equation: N ( u ) = −∇ · ( u (cid:48) ⊗ u (cid:48) ) ≡ f , (4)where the mean-pressure p is obtained such that ∇ · u = 0. The equation defining the mean-flow differs from thefixed point one ( N ( u ) = 0) by the presence of the right-hand side term f ≡ −∇ · ( u (cid:48) ⊗ u (cid:48) ). The goal here is tomodel this term by approximating the fluctuation u (cid:48) by singular vectors of the Resolvent operator. To do so, wesubtract the mean flow equation (4) from the instantaneous equations (1), leading to: ∂ t u (cid:48) + L u u (cid:48) = ∇ · ( u (cid:48) ⊗ u (cid:48) − u (cid:48) ⊗ u (cid:48) ) ≡ f (cid:48) , (5)where the operator L u u (cid:48) = u · ∇ u (cid:48) + u (cid:48) · ∇ u + ∇ p (cid:48) − Re − ∆ u (cid:48) is the linearized Navier-Stokes operator and f (cid:48) represents the nonlinear fluctuations, acting as a forcing term. We then try to find solutions of this equation inthe frequency domain by using the ansatz u (cid:48) = ˆ u e i ωt + c.c. and f (cid:48) = ˆ f e i ωt + c.c. , leading to:ˆ u ( ω ) = R ( u , ω )ˆ f ( ω ) , where R ( u , ω ) ≡ (i ωI + L u ) − (6)is the so-called Resolvent operator. This equation represents an input/output relation between the forcing termˆ f and its perturbation ˆ u . In order to isolate the most energetic dynamics (Beneddine et al. (2016)), we maximizethe gains: G ( ω ) = || ˆ u || || ˆ f || = ( R † ( u ; ω ) R ( u ; ω )ˆ f , ˆ f )(ˆ f , ˆ f ) (7)where R † ( u , ω ) is the adjoint Resolvent operator, given by the general definition ( R † a , b ) = ( a , R b ) , ∀ a , b . Here weconsider the energy inner product ( a , b ) = (cid:82) Ω a ∗ · b d x . This optimization problem may be solved by consideringeither of the following eigen-problems: R † ( u , ω ) R ( u ; ω )ˆ f k = µ k ( ω )ˆ f k , R ( u , ω ) R † ( u ; ω )ˆ u k = µ k ( ω )ˆ u k , (8)where µ k are the eigenvalues (ordered such that µ ≥ µ ≥ · · · ≥ u k and ˆ f k are the optimal responses andforcings, both normalized such that || ˆ f k || = || ˆ u k || = 1. These vectors consist in two orthonormal bases withrespect to the chosen inner-product and may be related through ˆ u k = µ − k R ( u ; ω )ˆ f k . A Fourier mode of theflowfield ˆ u can therefore be linked to the nonlinear forcing mode ˆ f through:ˆ u ( x , ω ) = + ∞ (cid:88) k =1 µ k ( ω ) β k ( ω )ˆ u k ( x , ω ) , (9)where β k ( ω ) = (cid:16) ˆ f k ( ω ) , ˆ f ( ω ) (cid:17) is the projection of the nonlinear forcing term onto the right singular-vector. Inshear-dominated flows, it is common that rank-1 approximations of the Resolvent operator may be sufficient torepresent the input-output dynamics of the flow, leading to the following approximation of the fluctuation field:ˆ u ( x , ω ) ≈ µ ( ω ) β ( ω )ˆ u ( x , ω ) , (10)At this point, the notations ˆ u ( ω ) and ˆ f ( ω ) refer to a given spectral decomposition of the flowfield: in the following,we will make a distinction between time-periodic (usually oscillators, for which β k ( ω ) is a Dirac-like distribution)and broadband-frequency (usually noise-amplifiers, for which β k ( ω ) is broadband) flows. .1 Time-periodic flows Many transitional flows present a periodic behaviour such as the flows around bluff bodies. In those cases thesteady-state solution is unstable and the instability converges toward a saturated periodic limit-cycle. Further-more, it is not uncommon that most of the energy of this periodic limit-cycle is concentrated at the fundamentalfrequency ω (see Turton et al. (2015)). For those flows, the harmonic averaging process:ˇ u ( ω ) = lim T → + ∞ T F T ( u (cid:48) )( ω ) , where F T ( u (cid:48) )( ω ) = (cid:90) T e − i ωt u (cid:48) ( t ) dt (11)converges for some discrete frequencies ω and the lowest value corresponds to the fundamental mode of theFourier series expansion (or, more generally, a Koopman mode of the flow, see Arbabi and Mezi´c (2017)). Mezi´c(2013) and Turton et al. (2015) showed that a Koopman mode was close to the marginal eigensolution of thelinearized Navier-Stokes operator at the frequency ω if the flow is harmonic (see Barkley, Gomes, and Henderson(2002); Mantic-Lugo et al. (2015); Sipp and Lebedev (2007)). Here we will approximate the Koopman modewith the Resolvent mode ˇ u ( ω ) ≈ A ˆ u ( ω ), which is justified by the fact that these structures coincide when theeigenvalue is almost purely imaginary. Furthermore, even though the gain µ ( ω ) is large, the nonlinear forcingcoefficient ( β ( ω )) is generally weak (interaction between the weak second harmonic at 2 ω and the fundamentalat ω ), yielding a finite amplitude A = µ ( ω ) β ( ω ). The fluctuation then reads: u (cid:48) ≈ A ˆ u e i ω t + c.c., (12)where the parameter A is now a single complex scalar. Based on this, we can approximate the Reynolds-stresstensor as f = − | A | (cid:60){ ˆ u ⊗ ˆ u ∗ } . Considering the mean flow equation (4), an approximation ˜ u of the mean-flow u may therefore be obtained through: N (˜ u ) = −| A | φ (ˆ u , ˆ u ∗ ) (13a) R (˜ u , ω ) R † (˜ u , ω )ˆ u = µ ˆ u , (13b)where ∗ stands for the complex-conjugation and φ ( u , u ) = ( u · ∇ ) u + ( u · ∇ ) u is the symmetrized convectiveoperator. We also remind that the dominant optimal response has been normalized such that || ˆ u || = 1. At thispoint, we remark that the system formed by the mean flow equation and the Resolvent analysis (13) is a closedset of nonlinear equations, if one knows the frequency of the nonlinear signal ω and its energy | A | (or amplitude | A | ). Arguably, the frequency of the flow may be determined from a time-resolved signal of the flow, for instancea hot-wire or a point-wise wall-pressure measurement. Furthermore, the amplitude of that signal (or any othermeasure of the flowfield) could be used to determine the amplitude of the resolvent mode | A | since, for everyvalue of | A | , one can solve for the nonlinear solution of (13) and therefore tune this parameter (if the relation ismonotonous) so that the solution matches the actual measure. The output of this procedure is the full flow-field u composed of the reconstructed mean-flow ˜u and the fluctuation u (cid:48) given in equation (12). The strength of theapproach lies in the simplicity of the optimization problem, since it only involves a single parameter contrarily tothe procedure described in Foures et al. (2014), where large-scale spatial fields had to be determined. For broadband flows, the previous approximation is no longer suited. Indeed, since the frequency spectrum isbroadband, we need to ”discretize” the fluctuation field in the frequency domain. Furthermore, the fluctuationfield u (cid:48) does not hold any Koopman mode (or, alternatively, Fourier-series representation), since the quantityˇ u ( ω ) converges to zero for all frequencies ω . For this reason, we consider the statistics of the unsteady flow,contained in the two-point ( x and x (cid:48) ), two-time (the configuration being homogeneous in time, a single time t isrequired) correlation tensor: T x , x (cid:48) ( t ) := E [ u (cid:48) ( x , τ ) u (cid:48) ( x (cid:48) , τ + t )] = lim T → + ∞ T (cid:90) T u (cid:48) ( x , τ ) u (cid:48) ( x (cid:48) , τ + t ) dτ. (14)Due to the absence of Koopman modes, this tensor tends to zero when t → ±∞ and is therefore square-integrable.For this reason, we may consider the Fourier-transform of this quantity to obtain the two-point spectral correlationtensor: ˆ T x , x (cid:48) ( ω ) := (cid:90) + ∞−∞ e − i ωt T x , x (cid:48) ( t ) dt = lim T → + ∞ T E [ F T ( u (cid:48) ( x ))( ω ) F T ( u (cid:48) ( x (cid:48) ))( ω ) ∗ ] (15)If x = x (cid:48) , T x , x ( t ) is called the auto-correlation function such that, for t = 0, T x , x (0) contains the Reynolds-stress tensor information at x . ˆ T x , x ( ω ) is the power-spectral density (PSD), which contains the distribution ofthe energy of the flow per frequency at x . Using the approximation by the Resolvent analysis (10), this tensorbecomes: ˆ T x , x (cid:48) ( ω ) ≈ µ E [ | β | ]ˆ u ( x , ω )ˆ u ∗ ( x (cid:48) , ω ) , (16) T x , x ( ω ) for a broadband flow (a). Its integral onfive ( N = 5) frequency intervals (represented as the shades of gray) are compared with the discrete-in-frequency energies coming from the model (b). where, since the SVD decompositon is deterministic (since it is based on the mean-flow quantity), the expectationoperator acts only on the expansion coefficients β k ( ω ). Coming back now to the cross-correlation tensor, applyingthe inverse Fourier Transform on this tensor, we have: T x , x (cid:48) ( t ) = 12 π (cid:90) + ∞−∞ ˆ T x , x (cid:48) ( ω ) e i ωt dω ≈ π (cid:90) + ∞ µ E [ | β | ]ˆ u ( x , ω )ˆ u ∗ ( x (cid:48) , ω ) dω + c.c. (17)This approximation is still not adequate for our purpose since it is not yet discrete in frequency. To do so, weconsider a ”grid” in the frequency space, { Ω j } j =1 ,N +1 , so that most of the energy of the fluctuation falls in theinterval (Ω , Ω N +1 ). Furthermore, in the following, we set Ω = 0. A further approximation of (17) is finallygiven by: T x , x (cid:48) ( t ) ≈ N (cid:88) j =1 | A j | ˆ u ( x , ω j )ˆ u ∗ ( x (cid:48) , ω j ) + c.c., (18)where we conveniently define | A j | = (2 π ) − µ ( ω j )(Ω j +1 − Ω j ) E [ | β ( ω j ) | ] and the frequencies ω j are definedto be the center point of each of the intervals in frequency ω j = (Ω j + Ω j +1 ) /
2. This approximation leads to afluctuation u (cid:48) of the following form: u (cid:48) = N (cid:88) j =1 A j ˆ u ( ω j ) e i ω j t + c.c. (19)which is nothing but an extension of the approximation employed for the time-periodic flow to consider multiplefrequencies. The approximation of the mean-flow ˜ u can be given now using equation (4), implying that: N (˜ u ) = − N (cid:88) j =1 | A j | φ (ˆ u ( ω j ) , ˆ u ∗ ( ω j )) (20a) R ( u , ω j ) R † ( u , ω j )ˆ u ( ω j ) = µ ( ω j )ˆ u ( ω j ) , (20b)where again all dominant optimal responses are unit-norm.Similarly as before, this system is closed if we provide the amplitudes | A j | of the Resolvent modes. Again,we will suppose that we have access to (few) point-wise time-resolved probes in the flow, located at some probingpoints { x mp } m =1 , ··· ,M , allowing us to have a good idea not only of the amplitude | A j | (and its phase) but also ofthe frequencies ω j . We need then to establish a comparison basis between the “discrete-frequency” model 20 andits “continuous-frequency” counterpart, given by the real flow. To do so, we integrate equation (16), evaluated atthe measurement points x mp = x = x (cid:48) , leading to: I x mp ,ω j ≡ (cid:90) Ω j +1 Ω j ˆ T x mp , x mp ( ω ) dω ≈ π | A j | | ˆ u j ( x mp , ω j ) | . (21)In Figure 1, we provide a schematic view of this integral: the PSD is divided in a few intervals (resulting nowin an energy quantity) and “lumped” at the discrete set of frequencies { ω j } .In terms of practical resolution of the system, the inputs are the quantities I x mp ,ω j and the scanning of thefull control space {| A j |} ∈ R N is no longer feasible as soon as N becomes large. For this reason, we will adopt anoptimization approach where we will penalise the differences between I x mp ,ω j and 2 π | A j | | ˆ u j ( x mp , ω j ) | for somemeasurement points and frequency bands. We remark that, in theory, only one point ( M = 1) is sufficient fordetermining all the amplitudes. However, in practise, that point could have weak energies for certain range offrequencies, leading often to poor predictions (as discussed in Beneddine et al. (2016)). For this reason, we willchoose specific frequency bands within the list of available measurement points. More details on the constructionof the cost functional will be given later. Once this cost-functional is chosen, its minimization will be carried outby a gradient-free method, here the COBYLA algorithm (see Conn, Scheinberg, and Toint (1997); Powell (2007)),from python/scipy’s library “minimize”. It is a simplex-based algorithm, which is appropriate if the number ofmodes N remains reasonably small. For larger N (or if sub-optimal response modes need to be considered) agradient is needed, which can be obtained from an adjoint-based approach. .2.1 Spectral analysis of broadband flows At this point, it is convenient to present the tools to perform a spectral analysis of broadband flows. Theresulting quantities will be compared with the mean-flow Resolvent modes, as done later to assess the efficiencyof the fluctuation reconstruction. This analysis is based on the decomposition of the spectral cross-correlationtensor ˆ T x , x (cid:48) ( ω ), which may hold several spatio-temporal coherent modes that can be extracted, from the mostenergetic to the least. To do so, we employ the Spectral Proper-Orthogonal Decomposition (SPOD, see Towne etal. (2018)): (cid:90) Ω ˆ T x , x (cid:48) ( ω ) ˆ ψ k ( ω, x (cid:48) ) d x (cid:48) = σ k ( ω ) ˆ ψ k ( ω, x ) , (22)where the eigenvalues σ k indicate the strength of the coherent (SPOD) mode ˆ ψ k . This tensor is computed withthe Welch method where the expectation operator E [ · ] is approximated by decomposing the signal of a single run(DNS or experimentally obtained) into N b overlapping and sufficiently long bins (of length T ), leading to:ˆ T x , x (cid:48) ( ω ) ≈ T N b N b (cid:88) k =1 F kT ( u (cid:48) ( x (cid:48) ))( ω ) F k, ∗ T ( u (cid:48) ( x ) ∗ )( ω ) ≡ ˆ U ˆ U ∗ ,T , (23)where F kT is the previously defined truncated Fourier transform, acting on the k -th bin, and ˆ U ( ω ) = [ · · · , F kT ( u (cid:48) )( ω ) , · · · ] / √ T N b is the (normalized) collection of Fourier-Transforms of all the bins. The superscript T denotes transposition. Ina discrete in space system, the eigensystem becomes then: U ( ω ) U ∗ ( ω ) B ˆ ψ k ( ω ) = σ k ( ω ) ˆ ψ k ( ω ) , (24)where B is a mass-matrix, containing the mesh-metrics, required by the space integral in (22). This eigensystem,if too large to be easily solved, can be transformed in a smaller one by considering the eigensystem for the variableˆ y k : ˆ U ∗ ( ω ) B ˆ U ( ω )ˆ y k ( ω ) = σ k ( ω )ˆ y k ( ω ) , (25)and the SPOD mode can be recovered from ˆ ψ k ( ω ) = σ − k ˆ U ( ω )ˆ y k ( ω ). Generally speaking, the dominant SPODmode should be comparable with the dominant optimal response mode if the Resolvent operator is close to rank1, that is σ (cid:29) σ (Beneddine et al., 2016). If the expansion coefficients β k ( ω ) are uncorrelated E [ β k β l ] = αδ k,l ,which is a strong requirement, then all SPOD modes should correspond to the optimal response modes (Towneet al., 2018). In this paragraph, we describe the numerical method to solve the system of nonlinear equations given by (13)or (20). One possible way is to iteratively solve the mean flow for a given dominant optimal response mode andthen solve for a new mode with the new mean flow and keep alternating up to the point where a fixed pointsolution is reached. This method was employed by Mantic-Lugo et al. (2015) and Mantic-Lugo and Gallaire(2016) on a similar model. However, it was observed that for relatively high values of | A | , the method requires alow under-relaxation factor, making the convergence slow. Here we propose an alternative method, which consistsin considering the set of equations 13 as a whole and applying a Newton method to it. To do so, we first rewritethe eigensystem determining the dominant optimal response with the variables ( λ = µ , ˆ u , ˆ a = λ − ˆ f ): λ (i ω B + L ˜ u )ˆ u = B ˆ a , ˆ u ∗ ,T B ˆ u = 1(i ω B + L ˜ u ) ∗ ˆ a = B ˆ u , (26)where ∗ stands for complex-conjugation, T stands for vector-transposition and B is the mass-matrix representingthe energy scalar-product. We remark here that, for the sake of clarity, we keep the ‘continuous’ notationintroduced in previous paragraphs, even if now those objects are discrete, since we are dealing with the numericalmethod. Together with those equations, we need to solve the mean-flow equation as well (for simplicity, we onlyshow the single-frequency case given by equation 13). Its linearization, together with the linearization of therewritten eigensystem leads to: L ˜ u | A | φ (ˆ u ∗ , ( · )) + | A | φ (ˆ u , ( · ) ∗ ) 0 0 λφ (( · ) , ˆ u ) λ (i ω B + L ˜ u ) (i ω B + L ˜ u )ˆ u − B u ∗ ,T B ( · ) 0 0 φ ∗ ,T (( · ) , ˆ a ) − B ω B + L ˜ u ) ∗ ,T δ ˜ u δ ˆ u δλδ ˆ a = − N ( ˜u , ˜u ) − | A | φ (ˆ u , ˆ u ∗ ) B ˆ a − λ (i ω B + L ˜ u )ˆ u (1 − ˆ u ∗ ,T B ˆ u ) / B ˆ u − (i ω B + L ˜ u ) ∗ ,T ˆ a . (27)This linear system is then solved with a GMRES iterative solver, preconditioned with the lower triangularmatrix (lower Gauss-Seidel), completing the description of the numerical method used. L .Table 1: Mesh dependency of steady and unsteady simulations for square cylinder at Re = 100 - chosenmesh for the forthcoming computations in boldface. Observed parameters are recirculation bubble’slength L , nonlinear frequency ω , the norm of the first harmonic energy || ˇ u ( ω ) || and dominant gain,evaluated at the fundamental frequency of the flow µ ( ω ).Steady Solution DNS L µ ( ω ) L ω || ˇ u ( ω ) || Mesh 1 43000 8.42 213775 2.43 0.913 2.131
Mesh 2
The two-dimensional flow around a squared-section cylinder is a typical example of oscillator flow. For Reynoldsnumber larger than Re ≈
50 (see Sohankar, Norberg, and Davidson (1998)), an unsteady behaviour naturally setsin, corresponding to a periodic limit cycle, characterized by a single frequency and its harmonics. Below, we willfirst ( § § A sketch of the computational domain is provided in Figure 2. The inflow boundary is located 15 diametersupstream of the center of the cylinder and a uniform velocity profile parallel to the lateral sides of the square isprescribed. The Reynolds number is set to Re = DU ∞ /ν = 100. The lateral boundaries are located 20 diametersaway from the center of the cylinder, on which we impose symmetry boundary conditions. The outflow boundarycondition is imposed 30 diameters downstream of the cylinder and reads ( pI − Re − ∇ u ) · n = . The referenceflow, that will be probed, is obtained by time-averaging the solution of a Direct Numerical Simulation (DNS). Thespatial discretization of both the DNS and the model relies on a Finite-Element Method (FEM) implemented inthe FreeFEM++ code (see Hecht (2012)). We use second-order Taylor-Hood elements ( P for velocity and P forthe pressure fields). The DNS solver employs a second-order semi-implicit temporal scheme for time-advancement.In Table 1, we show various quantities of interest to characterize the mean- and unsteady features of theflowfield: the recirculation bubble’s length L (for the steady solution and the time-averaged unsteady solutions),the fundamental frequency ω (obtained via a Discrete-Fourier Transform on a sufficiently long DNS run) andthe amplitude of the Fourier mode at the fundamental frequency in the unsteady simulation || ˇ u ( ω ) || . Resultsare provided for 3 meshes characterized by different grid densities. All of them are strongly refined around thesquare and in its wake and are coarsened in the free-stream region. We can see that mesh 2 provides a goodagreement with the finer mesh 3 for all observed quantities. Mesh 2 will therefore be chosen as the default meshin all following computations. In particular, the model will be computed on this mesh. Let us assume that we only know the time evolution of (for example) the cross-stream velocity component at asingle point in the wake of the cylinder, for example x m = ( x, y ) = (2 , u y at ( x, y ) = (2 , ω = 0 . | ˇ u y | ( ω ) = 0 .
21. Tuning of | A | in model (13) with this quantity (c), the optimal value being | A | ≈ . Re = 100. Stream-wise component of (a,b) mean-flow velocities (c,d)first-harmonic streamwise velocities A ˆ u and ˇ u ( ω ) and (e,f) mean force induced by the Reynolds stresstensor ˜ f = 2 | A | ∇ · (cid:60){ ˆ u ⊗ ˆ u ∗ } with | A | = 2 .
22. The black curves delimit the recirculation regions of thereconstructed mean flow. represents its harmonic-averages. In Figure 3 (b) we plot its modulus | ˇ u y | as function of frequency. This quantityindicates that the fundamental mode at ω strongly dominates the harmonics pω with p ≥
2. This suggests thata reconstruction with a single frequency may be sufficient. The fundamental frequency of the periodic limit cyclesis equal to ω = 0 .
91, while the amplitude is | ˇ u y ( x m , ω ) | = 0 . | A | in model (13), we may explore the control parameter space by com-puting the reconstructed solution for a range of | A | and therefore tuning | A | with the reference measurement datashown in figures 3(a,b). From a numerical point of view, this does not necessarily represent a high computationaleffort since one may proceed by progressively increasing the value of | A | and restart each new computation fromthe previous solution (for | A | = 0, the reconstructed mean-flow solution is the base-flow). Changing the value of | A | in sufficiently small steps, the reconstructed solution is obtained in a few Newton iterations. As | A | increases,the number of GMRES iterations to invert eq. (27) increases since the preconditioner becomes less accurate.Once the solutions of the model are obtained, we compare the prediction | A || ˆ u ( x m , ω ) | (shown with a solidline in figure 3(c)) with the measurement | ˇ u y | ( x m , ω ) (horizontal dashed-line). We remark that the predictedquantity is always increasing with | A | , meaning that there exists a unique value of | A | for which the predictedand measured values match. This procedure leads to a value of | A | = 2 .
22. This value slightly overestimatesthe actual DNS values, since || ˇ u ( ω ) || = 2 .
13. This is due to the fact that the model only takes into accountthe fundamental harmonic and not the higher order harmonics. The effect of the latter modes on the mean-flowdistortion is compensated by a slight overestimation of the amplitude of the fundamental Fourier mode. However,despite this small overestimation of the mode’s energy, the reconstructed solution of the model 13, representedin Figure 4 by the mean-flow (a), Fourier mode (b) and Reynolds-stress divergence (c), compares very well withthe reference counterparts (b,d,f). Furthermore, if we consider other flow quantities (such as the recirculationlength L , the mean-drag C D or the error related to the mean-flow estimation e = || ˜ u − u || ) to tune | A | , we cansee in Figure 5 that all of them yield a similar value | A | , between 2 .
22 and 2 .
25. This shows the robustness ofthe approach for this particular flow, which stands as an appealing alternative to the more classical mean-flowdata-assimilation approaches (Foures et al., 2014), especially if only few scalar measurements are available.At this point, we remind that a similar model was proposed by Mantic-Lugo et al. (2015) for the flow around acircular cylinder. Two major differences can be highlighted. First, instead of using optimal response modes, theyuse global modes for the approximation of the harmonics. Second, instead of considering external data to providea criterion for the choice of | A | , their model searched for a value of | A | such that the eigenmode was marginallystable. This condition is justified for flows behaving harmonically, where higher order harmonics are negligiblewith respect to the fundamental. Hence, their method is almost exclusive to harmonic ”oscillator” flows, where e = (cid:107) ˜ u − u (cid:107) . For each different measure, we may have a different(but close) value of the estimated | A | .Figure 6: Sketch of the physical domain and inflow boundary condition for the backward-facing step case.The length of the main recirculation bubble, the detachment and reattachment points at the top wall areindicated in the figure and defined to be L , x d and x r , respectively. unsteadiness is triggered by intrinsic dynamics. In the next section, we show the generality of our approach byconsidering the case of an ”amplifier flow”, here a backward-facing step flow where the fluctuation is triggered byexternal noise, amplified through linear convective instabilities (Marquet and Sipp (2010)). In this section, we consider the (more challenging) backward facing step flow configuration, which exhibits abroadband behaviour. We first briefly present the configuration ( § The configuration is the Backward-facing step described in Barkley et al. (2002); Herv´e et al. (2012). The Reynoldsnumber, based on the height of the step H and the maximum inlet velocity U ∞ of the Poiseuille flow imposed atthe inlet, is fixed at Re = U ∞ H/ν = 500. The inflow boundary condition is located at 5 H upstream of the stepand the outflow is located at 50 H downstream, where, again, the boundary condition ( pI − Re − ∇ u ) · n = isimposed. The other boundaries of the domain correspond all to solid walls on which a no-slip boundary conditionis imposed. A sketch of the computational domain is given in Figure 6.For the present flow conditions, there exists a globally-stable steady-state solution (see Barkley et al. (2002),Blackburn, Barkley, and Sherwin (2008)) which exhibits strong linear selective-frequency amplification mecha-nisms (Marquet and Sipp (2010)). It can be triggered by considering upstream stochastic noise, whose amplitudeis tuned so that the dynamics of the perturbations is nonlinear and generates a mean-flow deformation. In thepresent study, we pick that forcing term to be: f ( x , t ) = C (cid:18) φ ( x ) (cid:19) w ( t ) , (28)where C = 10 is the amplitude of the signal (chosen sufficiently high so that nonlinear effects are effective), w ( t )is a (pseudo) white-noise and φ ( x ) = e −| x − x G | / σ is a Gaussian function centered at x G = ( − . , .
25) of width σ = 0 . t = 0 . ∼ Similarly to what was done in the time-periodic case, we start by considering time-series of signals at a few locationsof the domain. This shows roughly how the energy of the fluctuations is distributed in the domain, especially L , x d , x r ), and fluctuation energy for, the unsteady solution A = (cid:113) T (cid:82) T | u (cid:48) | . Chosen mesh for the computations in boldfaceSteady Solution Unsteady Solution L x d x r L x d x r A Mesh 1 35000 10.875 8.688 17.508 7.120 5.781 9.912 1.07
Mesh 2 x p = (0 , . x p = (4 , . x p = (10 , .
2) and (d) x p = (15 , . as a function of frequency. In Figure 7, we show the PSD (computed using the Welch method, considering atime-series of length 5000 time units, subdivided in 99 bins with 50% overlap, with a sample time of ∆ t = 0 . ω ∈ (0 ,
2) (figures b,c,d). Furthermore, the amplitude of those spectrais much higher than in (a), indicating that the dynamics at these frequencies stems from convective amplificationmechanisms linked to the shear-layer generated at the step. This is in accordance with Marquet and Sipp (2010),who showed, thanks to a Resolvent analysis, that the peak of energy amplification is around ω ≈ .
5. For thisreason, for our model, we will choose Ω N +1 = 2 and a uniform frequency discretization of ∆ ω = Ω j +1 − Ω j , forwhich ω j = ∆ ω ( j + 1 / x mp to the onesmeasured. The associated error measure is given by: J ω j , x mp = (cid:32) π | A j | | ˆ u j ( x mp ) | I x mp ,ω j − (cid:33) . (29)We are now able to define our global cost functional, as a sum over at least N of those “elementary” contributions J ω j , x mp , possibly choosing different points x mp and different frequencies ω j . In what follows, we will choose always N of those quantities so as to be able to tune the N control parameters | A | j . For robustness reasons, we pickpoints x mp and frequencies j so that fluctuation energy is high in the frequency band linked to ω j . For this reason,we will ignore the signal at x p where the fluctuation level is weak for all frequencies. Furthermore, we can seethat the signal x p , contrarily to x , p contains a significant amount of energy within 1 ≤ ω < N +1 . For thisreason, we will “trust” this signal in that frequency range. For the range Ω = 0 ≤ ω <
1, we will use either ofthe two signals x , p . These choices result in the following 2 cost-functionals that will be minimized: J = (cid:88) ω j ≥ J ω j , x p + (cid:88) ω j < J ω j , x p , J = (cid:88) ω j ≥ J ω j , x p + (cid:88) ω j < J ω j , x p . (30)We will consider different numbers N of frequency bands and analyse the effect of this parameter on thereconstructed field. We will consider the cases N = 2 , ,
10. Also, to speed up the convergence for the cases N = 5 and 10, we first perform the optimization for the case N = 2, which is quick. Then, we initialize themean-flow for N = 5 with the solution obtained for N = 2 and the | A j | coefficients by computing the Resolventmodes (again based on the mean-flow determined with N = 2) and setting | A j | = (cid:113) I x mp ,ω j / π/ | ˆ u j ( x mp ) | , where x mp is the current measurement point. We repeat the same procedure to initialize the optimization for N = 10based on the mean-flow obtained for N = 5. In figure 8 we represent the convergence of the optimization process for the case N = 5 and for the two costfunctionals J and J described above. We can see that starting from the case N = 2 provides indeed a goodinitial guess, especially for J , since the optimization did not drift apart too much from the initialization and theoptimisation did not take more than 20 iterations to converge up to a precision of J n /J ≈ − . For the cost | A j | and of the cost functional J n /J for J = J (solid lines)and J = J (dashed lines). Case N = 5Table 3: Reconstructed parameters in function of number of frequencies considered N . L x d x r e = || ˜ u − u || A = (cid:113) T (cid:82) T (cid:82) Ω | u (cid:48) | DNS 7 .
07 5 .
62 9 .
99 – 1 . J J J J J J J J J J N = 2 7.42 9.39 5.91 7.57 9.31 13.34 0.10 0.70 1.20 0.81 N = 5 6.25 8.67 5.24 6.93 8.16 12.01 0.05 0.40 1.33 0.92 N = 10 6.02 7.78 5.35 6.77 7.80 10.40 0.09 0.20 1.41 1.11 functional J , the optimization needed slightly more iterations, mainly because of a poor estimate of the initialamplitude for ω = 1. Furthermore, we observe that the estimate of the amplitude for ω = 0 . ω = 0 . J . The rate of convergence of the cost functional J n /J is almostexponential for both cost functionals.In table 3, we have compared, for the 2 cost functionals and 3 numbers of frequencies ( N = 2 , , N ), we resolve better the frequency interval and we tend to add more energy intothe system. This energy can be measured by the quantity A = (cid:113) T (cid:82) T (cid:82) Ω | u (cid:48) | , homogeneous to an amplitude,which for the model we use can be recast as A = (cid:113)(cid:80) j | A j | . We can also see that, for the cost functional J ,this value gets closer to the exact one (computed from the DNS snapshots). This is not the case for the firstcost functional, which over-predicts this quantity, especially as N increases. We can thus see that, although themeasurement points were chosen so that their signals contain a significant and relevant amount of energy, theypredict slightly different scenarios. A possible explanation of this fact can be given in figure 9 where we plot theenergy of the resolvent modes and of the SPOD modes along the measurement line y = 0 . x , showing that, for example, at ω = 0 .
2, the SPOD mode (dashed line) is stronger thanthe resolvent mode (solid line) for x <
12 and the other way round for x > .
5, suggesting that if the signal from x p ( x = 10) is used (which is the case for J ), the energy for this frequency is probably going to be overpredicted,increasing A . On the other hand, if x p ( x = 15) is used ( J ) this energy will tend to be more accurately predicted.Looking now at the mean-flow quantities, such as the detachment / reattachment lengths L, x d , x r , we cansee that, by increasing N , and consequently A , we have almost always a shortening of those lengths, with anexception for x d with J , where it stayed at a constant value, near the reference one. This shows that, generally,the more fluctuation energy is provided to the system, the shorter the bubbles. However, those quantities aloneare not capable to infer, in a conclusive manner, how well the mean-flow is reconstructed. For this reason, weevaluate as well the mean-flow error e = || ˜ u − u || . We can see that, although the second cost functional J provides a better prediction of the fluctuation energy A , its mean-flow prediction is poorer than with J , forwhich the mean-flow error is of the order of e ∈ (0 . , .
1) (the minimal value being observed for N = 5), against e ∈ (0 . , .
7) for J . This interplay between a good prediction of the energy of the fluctuation field and a goodmean-flow prediction was already observed in the previous cylinder configuration (although much more mitigated)where the value of A minimizing the mean-flow error e was higher than the DNS one due to the absence in themodel of the second-order harmonic. A similar argument can therefore be made here.We now turn our attention to the comparison of the assimilated fields with the ones from the DNS. In figure 10,we compare the reconstructed mean-flow and Resolvent modes for J and N = 5 with the mean-flow and SPODmodes of the DNS. We can see that, although the quantitative features are not necessarily very accurate (length ofupper separation bubble, amplitude of fluctuations), the qualitative features of the mean-flow and SPOD modesare well captured. This shows that the present procedure is effective to reconstruct the overall features of themean- and unsteady characteristics of a flow. We remind also that only very few point-wise measurements havebeen considered for that purpose. y = 0 .
2, of the resolvent modes onthe mean flow (solid lines) and Fourier modes (dashed lines) at several frequencies (both modes arenormalized). (a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 10: Comparison between reconstructed (a,c,e,f,i,k) and reference (b,d,f,h,j,l) fields: mean-flow(a,b) and Fourier modes at ω = 0 . ω = 0 . ω = 1 . ω = 1 . ω = 1 . σ k ( ω ), (b) Resolvent gains µ k ( ω ) and (c) the (normalized) projectionbetween resolvent and SPOD leading modes as a function of frequency. A better understanding of the quality of the reconstruction can be gained by analysing the SPOD / resolventmodes and the projection coefficients between these modes. We have plotted in figure 11 (a) the two dominantSPOD gains as a function of frequency. We can see that a quite strong suboptimal branch exists, indicating thatthe dynamics is not really rank one. These suboptimal SPOD modes have not been considered in the presentarticle, but could in principle be. Figure 11 (b) shows that the flow exhibits a dominant resolvent mode, whichis at the origin of the dominant SPOD mode. Yet, the suboptimal resolvent gain is not very weak with respectto the dominant one, which may explain why the dynamics is not strictly rank-one in terms of SPOD gains. In11 (c) we can see that the optimal SPOD and resolvent modes are not fully aligned, a property which is at theheart of the present reconstruction procedure. Hence, in flowfields where this alignement is even better and thedynamics closer to rank 1 (see for example turbulent jets where the alignement reaches values of 0.9 instead of0.6 here, see Schmidt, Towne, Rigas, Colonius, and Br`es (2018)), the reconstruction could be much better. Notealso that a more complex model with additional suboptimal modes could be considered. Also, Pickering et al.(2020) showed that the alignement of the optimals (and the suboptimals) could be improved in turbulent flowsby considering an eddy-viscosity for the definition of the resolvent modes.
In this work, we presented a data-assimilation technique based on Resolvent modes around the mean-flow. Thetechnique was based on a model composed of a mean-flow equation coupled with resolvent modes correspondingto left singular vectors of the Resolvent operator. The model is not closed since the amplitudes of those modes areunknown and need to be tuned with the help of measurements, which are typically takes as point-ise time-resolvedprobes. This technique was applied on an oscillator flow, whose frequency content was essentially monochromatic.For this reason, the data-assimilation results showed that we are indeed capable of accurately recovering not onlythe mean-flow but also the fluctuation around it. This technique was also applied on a more challenging noise-amplifier backward-facing step flow where the dynamics is here driven by an external white-noise forcing, leadingto a broadband-frequency fluctuation. For this reason, a connection between this broadband fluctuation andour model, which is essentially peaked-frequency, needed to be established. This connection is based on theintegral of the Power-Spectral Density, which is a finite quantity either for the actual reference flow and for ourmodel. By enforcing our model to satisfy those quantities at certain point locations, where the velocity field wasmeasured, we were able to obtain a fairly good solution that reproduced this time both mean-flow quantities andthe fluctuation as well. One shortcoming of this approach was the dependency of the assimilated solution on thespatial location of the measurements. This comes from a mismatch between the SPOD modes and the Resolventones. Indeed, if one of those modes is stronger/weaker than the other in a specific region of the flow, we may beover/under predicting its energy, leading to erroneous results. Another shortcoming of this approach is that onlythe leading mode is taken into account in the model, leading to other possible models where several modes, for agiven frequency, are considered in the model. Apart from the higher computational complexity, one difficulty ofthis approach is the tendency of suboptimal Resolvent and SPOD modes to differ. For those reasons, generally,we believe that adding an eddy-viscosity could improve the results since it usually improve the comparability ofthose modes (see Morra et al. (2019); Pickering et al. (2020)).
References
A., G., D., H., & E., M. (2013). Inflow and initial conditions for direct numerical simulation based onadjoint data assimilation.
J. Comp. Phys. , , 480-497.Arbabi, H., & Mezi´c, I. (2017). Study of dynamics in post-transient flows using koopman mode decom-position. Physical Review Fluids , (12), 124402.Barkley, D. (2006). Linear analysis of the cylinder wake mean flow. Europhys. Lett. , , 750–756.Barkley, D., Gomes, M. G. M., & Henderson, R. D. (2002). Three-dimensional instability in flow over abackward-facing step. J. Fluid Mech. , , 167–190.13eneddine, S., Sipp, D., Arnault, A., Dandois, J., & Lesshafft, L. (2016). Conditions for validity of meanflow stability analysis. Journal of Fluid Mechanics , , 485–504.Beneddine, S., Yegavian, R., Sipp, D., & Leclaire, B. (2017). Unsteady flow dynamics reconstruction frommean flow and point sensors: an experimental study. Journal of Fluid Mechanics , , 174–201.Blackburn, H. M., Barkley, D., & Sherwin, S. J. (2008). Convective instability and transient growth inflow over a backward-facing step. J. Fluid Mech. , , 271–304.Conn, A. R., Scheinberg, K., & Toint, P. L. (1997). On the convergence of derivative-free methodsfor unconstrained optimization. Approximation theory and optimization: tributes to MJD Powell ,83–108.D’Adamo, J., Papadakis, N., Memin, E., & Artana, G. (2007). Variational assimilation of pod low-orderdynamical systems.
J. Turbul. , , 1-22.Duraisamy, K., Iaccarino, G., & Xiao, H. (2019). Turbulence modeling in the age of data. Annual Reviewof Fluid Mechanics , , 357–377.Evensen, G. (2009). Data assimilation: The ensemble kalman filter. (2nd ed.). Springer-Verlag.Ferrero, A., Iollo, A., & Larocca, F. (2020). Field inversion for data-augmented rans modelling inturbomachinery flows.
Computers & Fluids , , 104474.Foures, D. P. G., Dovetta, N., Sipp, D., & Schmid, P. J. (2014). A data-assimilation method forreynolds-averaged navierstokes-driven mean ow reconstruction. J. Fluid Mech. , , 404–431.Franceschini, L., Sipp, D., & Marquet, O. (2020). Mean-flow data assimilation based on minimal cor-rection of turbulence models: application to turbulent high-reynolds number backward-facing step. arXiv preprint arXiv:2006.12151 .G´omez, F., Blackburn, H., Rudman, M., Sharma, A., & McKeon, B. (2016). A reduced-order modelof three-dimensional unsteady flow in a cavity based on the resolvent operator. Journal of FluidMechanics , .He, C., Liu, Y., Gan, L., & Lesshafft, L. (2019). Data assimilation and resolvent analysis of turbulentflow behind a wall-proximity rib. Physics of Fluids , (2), 025118.Hecht, F. (2012). New development in freefem++. J. Num. Math. , , 251–266.Herv´e, A., Sipp, D., & Schmid, P. J. (2012). A physics-based approach to flow control using systemidentication. J Fluid Mech. , , 26–58.Hwang, Y., & Cossu, C. (2010). Amplification of coherent streaks in the turbulent couette flow: aninput–output analysis at low reynolds number. Journal of Fluid Mechanics , , 333–348.Le Dimet, F. X., & Talagrand, O. (1986). Variational algorithms for analysis and assimilation ofmeteorological observations: theoretical aspect. Tellus 38A , 97-110.Liu, C., Xiao, Q., & Wang, B. (2008). An ensemble-based four-dimensional variational data assimilationscheme. part i: Technical formulation and preliminary test.
Mon. Weather Rev. , , 3363–3373.Lorenc, A. C. (1986). Analysis methods for numerical weather prediction. Quart. J. R. Met. Soc. , ,1177–1194.Mantic-Lugo, V., Arratia, C., & Gallaire, F. (2015). A self-consistent model for the saturation dynamicsof the vortex shedding around the mean flow in the unstable cylinder wake. Phys. Fluids , .Mantic-Lugo, V., & Gallaire, F. (2016). Self-consistent model for the saturation mechanism of theresponse to harmonic forcing in the backward-facing step flow. J. Fluid Mech. , , 777-797.Marquet, O., & Sipp, D. (2010). Global sustained perturbations in a backward-facing step flow. IUTAMBookseries , .McKeon, B. J., & Sharma, A. S. (2010). A critical-layer framework for turbulent pipe flow. J. FluidMech. , , 336–382.Meldi, M. (2018). Augmented prediction of turbulent flows via sequential estimators. Flow, Turbulenceand Combustion , (2), 389–412.Mezi´c, I. (2013). Analysis of fluid flows via spectral properties of the koopman operator. Annual Reviewof Fluid Mechanics , , 357–378.Mons, V., Chassaing, J. C., Gomez, T., & Sagaut, P. (1986). Reconstruction of unsteady viscous flowsusing data assimilation schemes. J. Comp. Phys. , , 255-280.Morra, P., Semeraro, O., Henningson, D. S., & Cossu, C. (2019). On the relevance of reynolds stressesin resolvent analyses of turbulent wall-bounded flows. arXiv preprint arXiv:1901.04356 .Parish, E. J., & Duraisamy, K. (2017). A paradigm for data-driven predictive modeling using fieldinversion and machine learning. J. Comp. Phys. , , 758–774.Pickering, E., Rigas, G., Schmidt, O. T., Sipp, D., & Colonius, T. (2020). Optimal eddy viscosity forresolvent-based models of coherent structures in turbulent jets. arXiv preprint arXiv:2005.10964 .Pier, B. (2002). On the frequency selection of finite-amplitude vortex shedding in the cylinder wake.14 ournal of Fluid Mechanics , , 407–417.Powell, M. J. (2007). A view of algorithms for optimization without derivatives. Mathematics Today-Bulletin of the Institute of Mathematics and its Applications , (5), 170–174.Schmidt, O. T., Towne, A., Rigas, G., Colonius, T., & Br`es, G. A. (2018). Spectral analysis of jetturbulence. Journal of Fluid Mechanics , , 953–982.Singh, A. P., & Duraisamy, K. (2016). Using field inversion to quantify functional errors in turbulenceclosures. Phys. Fluids , .Sipp, D., & Lebedev, A. (2007). Global stability of base and mean flows: a general approach and itsapplications to cylinder and open cavity flows. J. Fluid Mech. , , 333–358.Sipp, D., & Schmid, P. J. (2016). Linear closed-loop control of fluid instabilities and noise-inducedperturbations: A review of approaches and tools. Applied Mechanics Reviews , (2).Sohankar, A., Norberg, C., & Davidson, L. (1998). Low-reynolds-number flow around a square cylinderat incidence: study of blockage, onset of vortex shedding and outlet boundary condition. Int.J. Numer. Methods Fluids , , 39–56.Symon, S., Dovetta, N., McKeon, B. J., Sipp, D., & Schmid, P. J. (2017). Data assimilation of meanvelocity from 2d piv measurements of flow over an idealized airfoil. Exp. Fluids , , 1–17.Symon, S., Sipp, D., & McKeon, B. J. (2019). A tale of two airfoils: resolvent-based modelling of anoscillator vs. an amplifier from an experimental mean. arXiv preprint arXiv:1904.10131 .Symon, S., Sipp, D., Schmid, P. J., & McKeon, B. J. (2020). Mean and unsteady flow reconstructionusing data-assimilation and resolvent analysis. AIAA Journal , (2), 575–588.Towne, A., Lozano-Dur´an, A., & Yang, X. (2020). Resolvent-based estimation of space–time flowstatistics. Journal of Fluid Mechanics , .Towne, A., Schmidt, O. T., & Colonius, T. (2018). Spectral proper orthogonal decomposition and itsrelationship to dynamic mode decomposition and resolvent analysis. Journal of Fluid Mechanics , , 821–867.Turton, S. E., Tuckerman, L. S., & Barkley, D. (2015). Prediction of frequencies in thermosolutalconvection from mean flows. Phys. Rev. E. ,
91 (4), 04300991 (4), 043009