domain: spi phys

text:
Colin Leclercq

Fabrice Demourant

Charles Poussot-Vassal

Denis Sipp

Linear iterative method for closed-loop control of quasiperiodic flows

 to study receptivity mechanisms in turbulent flows. Using this linear time-invariant model of the nonlinear flow, modern control methods such as structured H ∞ -synthesis can be used to design a controller. The approach is successful in damping self-sustained oscillations associated to specific eigenmodes of the mean flow spectrum. Despite excellent performance, the linear controller is however unable to completely suppress flow oscillations, and the controlled flow is effectively attracted towards a new dynamical equilibrium. This new attractor is characterized by a different mean flow, which can in turn be used to design a second controller. The method can then be iterated on subsequent mean flows, until the coupled system eventually converges to the base flow. An intuitive parallel can be drawn with Newton's iteration: at each step, a linearized model of the flow response to a perturbation of the input is sought, and a new linear controller is designed, aiming at further reducing the fluctuations. The method is illustrated on the well-known case of two-dimensional incompressible open-cavity flow at Re = 7500, where the fully developed flow is initially quasiperiodic (two-torus state). The base flow is reached after five iterations. The present work demonstrates that nonlinear control problems may be solved without resorting to nonlinear reduced-order models. It also shows that physically relevant linear models can be systematically derived for nonlinear flows, without resorting to black-box identification from input-output data; the key ingredient being frequency-domain models based on the linearized Navier-Stokes equations about the mean flow. Applicability to amplifier flows and turbulent dynamics has, however, yet to be investigated.

Introduction

Despite significant progress over the past decades, taming flow instabilities and turbulence remains an area of fundamental investigation [START_REF] Kim | A linear systems approach to flow control[END_REF][START_REF] Brunton | Closed-loop turbulence control: progress and challenges[END_REF]Sipp & Schmid 2016). The potential impact of flow control is however tremendous: drag reduction, lift increase, mixing enhancement, noise reduction, lean combustion to name a few industrial applications. Control strategies fall into two broad categories: active or passive, depending whether energy is expended to drive an actuator or not. Actuation may be used in an open-loop framework, in order to drive the system away from its natural state by exciting its nonlinearities. For this strategy to be effective, a large power input may be required. On the other hand, closing the loop has a direct impact on the intrinsic dynamics of a linear system, hence the actuation amplitude may be significantly reduced [START_REF] Cattafesta | Active control of flow-induced cavity oscillations[END_REF]. However, uncertainties in the form of noise or modelling errors may cause instabilities or performance loss of the closed-loop, which in turn render the approach ineffective.

1.1. On the relevance of linear methods for closed-loop control of nonlinear flows In the particular case of fluid flows, modelling issues arise from the dimensionality of the system and the strong nonlinearity of the governing equations. Dimensionality reduction is often achieved by projecting the dynamics onto a set of well-chosen modal structures, and well-established techniques exist to systematically find such bases, e.g. global modes, proper orthogonal decomposition (POD) modes, balanced modes (BPOD), dynamic modes (DMD), etc. [START_REF] Bagheri | Input-output analysis, model reduction and control of the flat-plate boundary layer[END_REF]; [START_REF] Semeraro | Feedback control of three-dimensional optimal disturbances using reduced-order models[END_REF]; [START_REF] Barbagallo | Closed-loop control of an open cavity flow using reduced-order models[END_REF]; [START_REF] Rowley | Model reduction for flow analysis and control[END_REF] for a comprehensive review). Alternatively, if system identification is used to derive a reduced-order model (ROM) directly from inputoutput data [START_REF] Cattafesta | Development of an adaptive weapons-bay suppression system[END_REF][START_REF] Kegerise | Adaptive identification and control of flow-induced cavity oscillations[END_REF][START_REF] Illingworth | Feedback control of flow resonances using balanced reduced-order models[END_REF][START_REF] Hervé | A physics-based approach to flow control using system identification[END_REF][START_REF] Dahan | Feedback control for form-drag reduction on a bluff body with a blunt trailing edge[END_REF][START_REF] Dalla Longa | Reducing the pressure drag of a d-shaped bluff body using linear feedback control[END_REF], the resulting models will typically be of manageable dimension. Suitable ROMs for closed-loop control may also be obtained by gradually increasing the spatial discretization from a coarse mesh [START_REF] Jones | Modelling for robust feedback control of fluid flows[END_REF].

Nonlinearity, on the other hand, remains a major issue as the standard control tool set relies extensively on the hypothesis of a linear time-invariant (LTI) plant. Early studies therefore focussed on the linear dynamics of small perturbations around a base flow, i.e. a fixed point of the Navier-Stokes equations. Techniques were successfully developed to suppress intrinsic instabilities or prevent extrinsic noise amplification: oscillator flows may be controlled by damping their unstable modes using a feedback controller, while amplifier flows may be controlled by rejecting incoming perturbations using an upstream sensor (see Sipp & Schmid (2016) for a review). In some cases, the designed controllers may also be robust to nonlinearities; for instance if the Reynolds number is weakly supercritical, in the case of an unstable base flow [START_REF] Flinois | Feedback control of unstable flows: a direct modelling approach using the Eigensystem Realisation Algorithm[END_REF], or if incoming perturbations are large, in the case of a stable base flow [START_REF] Hervé | A physics-based approach to flow control using system identification[END_REF]. However, for strong enough nonlinearity, linearization about the base flow becomes irrelevant and control fails (Schmid & Sipp 2016).

In the case of self-sustained oscillations, which are of interest in this paper, there are indeed (at least) two distinct dynamical equilibria: the steady fixed point, and the unsteady attractor. In the uncontrolled flow, the fixed point is linearly unstable and the attractor lies at a finite distance from it in phase space. In this situation, the goal of a controller is not just to locally stabilize the base flow but to ensure that any initial condition on the attractor is within the basin of attraction of the base flow. If the attractor is not too far away from the base flow, then a robust LTI controller based on some mean LTI model may suffice [START_REF] Li | Feedback control of combustion instabilities from within limit cycle oscillations using H∞-infinity loop-shaping and the ν-gap metric[END_REF]. But again, for strong enough nonlinearity, modelling the plant as a single LTI operator may be too crude an approximation, and even robust control methods based on it may fail.

Significant efforts were therefore directed towards nonlinear reduced-order modelling. Physics-based nonlinear ROMs usually rely on Galerkin projection of the full Navier-Stokes equations (as opposed to the linearized version around the fixed point), most often using POD modes. However, such models typically require fine tuning and calibration because of instabilities stemming from the structure of the expansion and its truncation [START_REF] Noack | Galerkin method for nonlinear dynamics[END_REF]. While designing nonlinear ROMs is a challenge of fundamental interest by itself, it is however not immediately clear that their enhanced fidelity can be made profitable for closed-loop control applications. Indeed, optimal control based on such ROMs represent an enormous computational load [START_REF] Bergmann | Optimal rotary control of the cylinder wake using proper orthogonal decomposition reduced-order model[END_REF] for realistic application in experiments. Besides, adjoint-looping only yields a specific control command optimized over a given time horizon for a given initial condition of the system.

An alternative to nonlinear Galerkin models is nonlinear identification from inputoutput data (2 nd order Volterra and bilinear filters, [START_REF] Pillarisetti | Adaptive identification of fluid-dynamic systems[END_REF]; NARX, [START_REF] Dandois | Narx modelling of unsteady separation control[END_REF]; NARMAX, [START_REF] Semeraro | Qualitative dynamics of wave packets in turbulent jets[END_REF], cluster-based ROM, [START_REF] Kaiser | Cluster-based reduced-order modelling of a mixing layer[END_REF]; SINDy, [START_REF] Loiseau | Constrained sparse galerkin regression[END_REF], etc.) but the resulting models are usually not used for feedback control. Artificial neural networks have been trained to learn complex input-output relations [START_REF] Lee | Application of neural networks to turbulence control for drag reduction[END_REF][START_REF] Efe | Control of subsonic cavity flows by neural networks-analytical models and experimental validation[END_REF], but a controller design simply based on inversion of the neural net does not bring any improvement compared to existing LTI controllers. Another novel approach to tackle nonlinearity is machine-learning control [START_REF] Duriez | Closed-loop control of experimental shear flows using machine learning[END_REF][START_REF] Gautier | Closed-loop separation control using machine learning[END_REF]: a model-free approach which seeks a nonlinear control law through trial-and-error. Although promising, the approach is not amenable to physical analysis and interpretation, which may hinder further improvements in case of failure.

Systematic methods for deriving nonlinear control laws, i.e. relations from outputs to inputs applicable at any time, although more useful, are in fact rarely implemented in practice. It has only been done so far with analytically-tractable 'mean-field' ROMs of order 4 at most, for flows dominated by at most two frequencies, i.e. a single one for the unforced flow plus one for the actuator [START_REF] King | Nonlinear Flow Control Based on a Low Dimensional Model of Fluid Flow[END_REF][START_REF] Luchtenburg | Turbulence control based on reduced-order models and nonlinear control design[END_REF][START_REF] Aleksić | Robust nonlinear control versus linear model predictive control of a bluff body wake[END_REF][START_REF] Aleksić-Roeßner | On the need of nonlinear control for efficient model-based wake stabilization[END_REF]. Furthermore, although superior in some cases [START_REF] Aleksić-Roeßner | On the need of nonlinear control for efficient model-based wake stabilization[END_REF], nonlinear controllers do not systematically show obvious improvements compared to their linear counterparts. More complex nonlinear ROMs are in fact generally linearized about their fixed point, in order to design traditional LTI controllers [START_REF] Samimy | Feedback control of subsonic cavity flows using reduced-order models[END_REF][START_REF] Nagarajan | Development and application of a reduced order model for the control of self-sustained instabilities in cavity flows[END_REF]. While this procedure removes a lot of information from the initial model, the physical meaning of the fixed-point is unclear, as it is not directly related to either the mean flow nor the base flow of the true system.

Finally, passivity-based methods are suited to finite-amplitude perturbations of the base flow and, therefore, hold great promise for the problem of nonlinear flow control [START_REF] Sharma | Model reduction of turbulent fluid flows using the supply rate[END_REF][START_REF] Sharma | Relaminarisation of Reτ = 100 channel flow with globally stabilising linear feedback control[END_REF][START_REF] Heins | Passivity-based output-feedback control of turbulent channel flow[END_REF]. The key idea is to see the fluid flow as a Lur'e system [START_REF] Khalil | Nonlinear systems[END_REF], i.e. the nonlinearity of the governing equations is a feedback forcing to the linearised dynamics about the base flow (more details will be given in §2.1 and 2.2). Since nonlinearity in the Navier-Stokes equations arises from conservative processes, it cannot contribute to the energy growth of external perturbations; i.e. it is passive. Therefore, one is only concerned with energy growth arising from the linear processes. Methods are available to design LTI controllers preventing linear energy growth altogether [START_REF] Sharma | Relaminarisation of Reτ = 100 channel flow with globally stabilising linear feedback control[END_REF] or, if impossible, at least minimising it [START_REF] Heins | Passivity-based output-feedback control of turbulent channel flow[END_REF], but these methods have only been applied to linearly stable flows so far.

This brief overview of available nonlinear flow control methods explains why, to date, the vast majority of successful studies rely on LTI models and controllers, even though the underlying plant is nonlinear.

1.2. A physics-based linear framework relevant to nonlinear flows:

mean-flow perturbation analysis

In the following, we distinguish the base flow, defined as a fixed-point of the Navier-Stokes equations, and the mean flow, which is the temporal average of a fully-developped flow.

Modal analysis

Since [START_REF] Hammond | Global dynamics of symmetric and asymmetric wakes[END_REF], many authors have considered mean-flow stability analysis as a means to retrieve both the nonlinear oscillation frequency of globally unstable base flows and the coherent structure oscillating at that fundamental frequency. Most studies were concerned with wake flows [START_REF] Hammond | Global dynamics of symmetric and asymmetric wakes[END_REF][START_REF] Pier | On the frequency selection of finite-amplitude vortex shedding in the cylinder wake[END_REF][START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF][START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF][START_REF] Mittal | Global linear stability analysis of time-averaged flows[END_REF][START_REF] Leontini | A numerical study of global frequency selection in the time-mean wake of a circular cylinder[END_REF][START_REF] Meliga | Sensitivity of 2-d turbulent flow past a d-shaped cylinder using global stability[END_REF][START_REF] Camarri | Stability analysis of experimental flow fields behind a porous cylinder for the investigation of the large-scale wake vortices[END_REF]Mettot et al. 2014b;[START_REF] Mantič-Lugo | Self-consistent model for the saturation mechanism of the response to harmonic forcing in the backward-facing step flow[END_REF][START_REF] Carini | Global stability and control of the confined turbulent flow past a thick flat plate[END_REF], where mean-flow distortion leads to a significant reduction of the length of the recirculation region behind the solid body, compared to the base flow. Other configurations have been considered as well, among which open-cavity flow [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF][START_REF] Meliga | Harmonics generation and the mechanics of saturation in flow over an open cavity: a second-order self-consistent description[END_REF], gas turbine fuel injectors [START_REF] Juniper | Absolute and convective instability in gas turbine fuel injectors[END_REF]) and thermosolutal convection [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF]. We note that in resonant flows, coherent structures are rigorously described by Koopman modes, and the oscillation frequencies correspond to Koopman eigenfrequencies [START_REF] Arbabi | Study of dynamics in post-transient flows using Koopman mode decomposition[END_REF].

Resolvent analysis

For globally stable flows, the temporal instability framework does not properly account for the noise amplifier dynamics. In parallel flows, spatial stability analysis may be used to predict the response to localized harmonic forcing. In non-parallel flows, the resolvent operator may be used instead, to predict the linear response of the flow to any spatially-distributed forcing. Connections between local and global analyses obviously exist for weakly non-parallel flows. Global spatial modes obtained from solving parabolized stability equations may approximate optimal response modes obtained by singular value decomposition (SVD) of the resolvent operator [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]. But the resolvent operator is more than a mere extension of spatial stability analysis to the global framework. It fully describes the input-output behaviour of the flow. In particular, SVD of the resolvent operator not only yields a set of optimal response modes approximated by spatial instability modes, it also provides the associated set of optimal forcing modes and optimal gains giving insight into the receptivity mechanisms at play in the flow.

When high-amplitude noise excites the nonlinearities in amplifier flows, coherent motions seem to extract energy from the mean flow, rather than the base flow [START_REF] Butler | Optimal perturbations and streak spacing in wallbounded turbulent shear flow[END_REF][START_REF] Chomaz | Global instabilities in spatially developing flows: non-normality and nonlinearity[END_REF][START_REF] Del Álamo | Linear energy amplification in turbulent channels[END_REF][START_REF] Cossu | Optimal transient growth and very large-scale structures in turbulent boundary layers[END_REF][START_REF] Pujals | A note on optimal transient growth in turbulent channel flows[END_REF]. And just like mean-flow temporal modes capture the dominant coherent structures of resonator flows, mean-flow spatial modes or optimal response modes (often called resolvent modes) seem to capture the dominant coherent structures of amplifier flows. In that context, the coherent structures may be rigorously defined as spectral POD modes [START_REF] Gudmundsson | Instability wave models for the near-field fluctuations of turbulent jets[END_REF][START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF][START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF], i.e. the set of orthogonal modes optimally correlated in both space and time [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF].

Resolvent analysis has been extensively used for investigating noise sources in turbulent jets [START_REF] Garnaud | The preferred mode of incompressible jets: linear frequency response analysis[END_REF][START_REF] Jeun | Input-output analysis of high-speed axisymmetric isothermal jet noise[END_REF][START_REF] Semeraro | Modeling of coherent structures in a turbulent jet as global linear instability wavepackets: Theory and experiment[END_REF][START_REF] Schmidt | Wavepackets and trapped acoustic modes in a turbulent jet: coherent structure eduction and global stability[END_REF][START_REF] Tissot | Sensitivity of wavepackets in jets to nonlinear effects: the role of the critical layer[END_REF]. The method, initially introduced in the context of pipe flow [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF]Sharma & McKeon 2013;[START_REF] Luhar | Opposition control within the resolvent analysis framework[END_REF], has also been applied to channel flow [START_REF] Moarref | A low-order decomposition of turbulent channel flow via resolvent analysis and convex optimization[END_REF][START_REF] Luhar | A framework for studying the effect of compliant surfaces on wall turbulence[END_REF][START_REF] Nakashima | Assessment of suboptimal control for turbulent skin friction reduction via resolvent analysis[END_REF], backwardfacing step flow [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]) and shock-wave boundary-layer over a bump (Sartor et al. 2015a) among other flows. Resolvent analysis of convectively unstable flows has also been used in order to identify the physical mechanisms underlying existing passive and open-loop control strategies for drag reduction. [START_REF] Jacobi | Dynamic roughness perturbation of a turbulent boundary layer[END_REF] and [START_REF] Mckeon | Experimental manipulation of wall turbulence: a systems approach[END_REF] considered the effect of dynamic roughness actuation on a zero-pressure gradient turbulent boundary layer in the resolvent framework. [START_REF] Luhar | Opposition control within the resolvent analysis framework[END_REF] considered the case of opposition control in pipe flow while [START_REF] Luhar | A framework for studying the effect of compliant surfaces on wall turbulence[END_REF] and [START_REF] Nakashima | Assessment of suboptimal control for turbulent skin friction reduction via resolvent analysis[END_REF] respectively studied the effect of compliant walls and suboptimal control by blowing and suction in channel flow.

Interestingly, the resolvent operator about the mean flow also appears to be relevant to globally unstable flows. Resolvent analysis has been applied to lid-driven cavity flow [START_REF] Gómez | A reduced-order model of three-dimensional unsteady flow in a cavity based on the resolvent operator[END_REF], open-cavity flow [START_REF] Liu | Resolvent analysis of compressible flow over a long rectangular cavity[END_REF], the buffeting transonic flow over an airfoil (Sartor et al. 2015b) (although in that case the perturbed flow is the fixedpoint of a RANS model, not the time-average of an unsteady RANS simulation) or the flow past a square cylinder [START_REF] Gómez | Data-driven approach to design of passive flow control strategies[END_REF]. In the latter paper, the authors analyzed the structure of the optimal forcing modes in order to design passive control devices able to optimally disrupt the self-sustaining process of wake oscillation. More recently, [START_REF] Liu | Resolvent analysis of compressible flow over a long rectangular cavity[END_REF] also used resolvent analysis to guide the design of an open-loop control strategy in open cavity flow. In essence, the resolvent operator about the mean flow seems to characterize the input-output behaviour of many fully-developped flows, regarless of their amplifier or resonator character. But despite compelling evidence of the relevance of this operator to control problems, this tool does not seem to have been included in the framework of closed-loop control yet [START_REF] Mckeon | Experimental manipulation of wall turbulence: a systems approach[END_REF].

Using the resolvent operator to control open-cavity flow in a closed-loop framework

In the present paper, we propose to control the fully-developed flow over a globally unstable open cavity, using the resolvent operator about the mean flow and a feedbackloop setup (see figure 1 and§3.1- §3.2 for a comprehensive description of the model). The configuration under study has a unique fixed point, the base flow, which is unstable at the operating Reynolds number of Re = 7500. As we shall see in this paper, the system naturally evolves towards a quasiperiodic regime with two fundamental incommensurate frequencies (2-torus). Our aim is to bring the system back to the base flow, using sole knowledge of the mean flow, a single probe, a single actuator and linear methods only. The resolvent operator about the mean flow will be used to model the transfer function of the fully-developed flow, while modal analysis will be used to interpret the results and guide the design of a robust linear controller, in the framework of structured H ∞ synthesis [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF]. As the controlled flow latches onto a new dynamical equilibrium, which is not the expected base flow, we will iterate the method on the new mean flow, and so on with subsequent mean states, until convergence. The choice of open-cavity flow is motivated by its frequent occurrence in the literature, as a convenient test-bed for validating new ideas; the cavity flow model used in this paper, in particular, has been used in many past studies (Rowley et al. 2006;[START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF][START_REF] Barbagallo | Closed-loop control of an open cavity flow using reduced-order models[END_REF][START_REF] Barbagallo | Input-output measures for model reduction and closed-loop control: application to global modes[END_REF][START_REF] Dergham | Model reduction for fluids using frequential snapshots[END_REF][START_REF] Illingworth | Feedback control of cavity flow oscillations using simple linear models[END_REF][START_REF] Poussot-Vassal | Parametric reduced order dynamical model construction of a fluid flow control problem[END_REF]Schmid & Sipp 2016;Sipp & Schmid 2016). Open-cavity flow also has many pratical applications, in particular in the aerospace field: weapon bays, wheel wells and fuselage openings for telescopes and sensors, etc. (Rowley & Williams 2006;[START_REF] Cattafesta | Active control of flow-induced cavity oscillations[END_REF].

The plan of the paper is as follows. In section §2, we introduce the theoretical framework for our iterative control strategy. In section §3, we present the cavity flow setup and describe the uncontrolled dynamics. The design methodology for the first controller is then explained in §4. The controlled flow is analyzed in §5 and a short discussion section follows in §6. Finally, the main results are recalled in the conclusion section §7.

Theoretical framework

In this section, we cast the governing laws of motion in an input-output form convenient for feedback control. We then introduce a novel iterative method for stabilizing linearly unstable flows with a sequence of linear models and controllers. The methodology for designing robust controllers is explained later in §4.

The Navier-Stokes equations as a Lur'e system

We start by considering the nondimensional, incompressible, Navier-Stokes equations governing the velocity and pressure fields u and p

∂ t u + (u • ∇)u = -∇p + Re -1 ∇ 2 u, ∇ • u = 0, (2.1)
without external forcing. The Reynolds number Re := U L/ν is based on characteristic velocity and length scales U and L of the problem, and on the kinematic viscosity ν of the fluid. Equations (2.1) are, of course, supplemented by appropriate boundary conditions and an initial condition. Next, we decompose the unknown field q := (u, p) T ( T denoting the conjugate transpose) into an arbitrary steady incompressible reference field Q := (U, P ) T satisfying the boundary conditions and a perturbation q ′ := (u ′ , p ′ ) T with zero boundary conditions, such that q = Q + q ′ . The unsteady perturbation q ′ is governed by the following time-evolution equation

B∂ t q ′ = [A Q + Ψ Q ]q ′ (2.2)
where

Bq ′ := ∂ t u ′ 0 , (2.3) A Q q ′ := -(U • ∇)u ′ -(u ′ • ∇)U + Re -1 ∇ 2 u ′ -∇p ′ ∇ • u ′ , (2.4) Ψ Q q ′ := -(U • ∇)U -∇P + Re -1 ∇ 2 U -(u ′ • ∇)u ′ 0 .
(2.5)

The Jacobian A Q and the operator Ψ Q respectively account for the linear and nonlinear parts of the evolution operator associated with the flow. Since Q is assumed to be time-

ψ Q q ′ t=0 q ′ R Q (s) + B Figure 2.
The Navier-Stokes equations as a Lur'e system. independent, the dynamical system is autonomous; A Q is a linear time-invariant (LTI) operator and Ψ Q is a static nonlinearity. Taking the Laplace transform of (2.2) leads to the relation

q ′ (s) = R Q (s)[Ψ Q q ′ (s) + Bq ′ t=0 ], (2.6)
where s is the Laplace variable, q ′ t=0 is the initial condition on the perturbation and

R Q (s) := (sB -A Q ) -1 (2.7)
is called the resolvent operator about the reference field Q. Such interconnection of an LTI operator R Q with a static nonlinearity Ψ Q , illustrated in figure 2, is called a Lur'e system in nonlinear control theory [START_REF] Khalil | Nonlinear systems[END_REF][START_REF] Sharma | Model reduction of turbulent fluid flows using the supply rate[END_REF][START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF]. The resolvent operator is defined for any s which is not an eigenvalue of the Jacobian operator, i.e. a solution of (sB

-A Q )q ′ = 0 (2.8)
with non-zero q ′ . Interestingly, the output q ′ remains bounded even if A Q has unstable eigenvalues, thanks to the presence of the nonlinear feedback operator. Moreover, we stress that the resolvent operator itself is well-defined along the imaginary axis as long as A Q has no marginally unstable eigenvalues. The choice of the arbitrary reference field Q is discussed in the next paragraph.

Open-loop and closed-loop transfer functions of fully-developped flows

In panels (a) and (c) of figure 3, we now consider the effect of external forcing due to an actuator and modelled by a time-independent spatial field B multiplied by a timevarying scalar amplitude n switched on at t = 0. We also introduce a sensor characterized by the time-independent pair (C, Y ) providing the local measurement

y(t) := C, q ′ (t) + Y,
(2.9)

where

q 1 , q 2 := D q * 1 • q 2 dD (2.10)
defines an inner product over the domain D and * denotes complex conjugation. The output y is linear with respect to the fluctuation q ′ and the constant Y represents an arbitrary offset. Panel (a) represents the case of open-loop forcing, while panel (c) represents a feedback-loop configuration with a total actuation amplitude u + n, where

u := Ky (2.11)
is due to an LTI controller of transfer function K(s). With external forcing, the output q ′ is respectively given by in the open-loop case, and

q ′ = R Q [Ψ Q q ′ + Bq ′ t=0 + Bn] (2.12) y ψ Q q ′ R Q (s) + n B + K(s) u + ψ Q q ′ t=0 q ′ R Q (s) + n B + (b) (c) (a) w + n w y + n K(s) + (d) G Q (s) G Q (s) B q ′ t=0 B y C, • 1/s + Y y C, •
q ′ = (I -R Q BK C, • ) -1 R Q Ψ Q q ′ + Bq ′ t=0 + Bn + BK Y s (2.13)
in the closed-loop case. The latter expression can be rewritten as

q ′ = R K Q Ψ Q q ′ + Bq ′ t=0 + Bn + BK Y s , (2.14)
where

R K Q (s) := (sB -A K Q (s)) -1 (2.15)
is the resolvent operator associated with the Jacobian operator

A K Q (s) := A Q + K(s)B C, . . (2.16)
of the coupled-system fluid/controller. Because of the nonlinear feedback term, the concept of a transfer function from n to q ′ is ill-defined. However, using the resolvent operators introduced above and lumping all the terms independent of n in an output disturbance term w, we can recast the governing equations (2.12) and (2.14) in a form that resemble input-output relations:

q ′ = R Q Bn + w and q ′ = R K Q Bn + w (2.17)
in the open-loop and closed-loop cases respectively. Similar to the unforced case, the output q ′ is bounded despite the potential presence of unstable eigenvalues in A Q or A K Q , thanks to the nonlinear feedback term (included in w) ensuring saturation. An important point to note at this point is that the definition of the Lur'e system associated with the flow depends on the arbitrary choice of a steady reference field Q. The passivity theorem states that if the two interconnected blocks R K Q and Ψ Q are passive, i.e. can only store or dissipate energy, then the fixed-point of the closed-loop is stable with respect to perturbations of any amplitude. For Q = q b , the static nonlinearity

Ψ q b q ′ = -(u ′ • ∇)u ′ 0 (2.18)
happens to be conservative, hence passive. Therefore, a possible strategy to reach the base flow from the regime of fully-developped oscillations consists in designing an LTI feedback controller making the closed-loop linear dynamics about the base flow passive. That means designing a controller K suppressing both modal growth and transient growth arising from the non-normality of the Jacobian operator A q b . This strategy has been pursued in recent works by [START_REF] Sharma | Relaminarisation of Reτ = 100 channel flow with globally stabilising linear feedback control[END_REF] and [START_REF] Heins | Passivity-based output-feedback control of turbulent channel flow[END_REF], for globally stable flows (i.e. no modal growth). However, as discussed in the introduction, the resolvent operator about the mean flow appears to be a better choice for describing the behaviour of unforced flows in the fully-developped regime. [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF] therefore proposed to consider the mean flow q as an appropriate alternative choice of reference field Q in (2.6), in order to define the unforced Navier-Stokes as a Lur'e system. In that case, the nonlinear output

Ψ q q ′ = (u ′ • ∇)u ′ -(u ′ • ∇)u ′ 0 (2.19)
was interpreted as an endogenous forcing term. We decide to follow the same approach here and also choose

Q := q (2.20)
to model the flow receptivity to the exogenous forcing field B.

In the case of unforced flows, the mean flow is unambiguously defined as the temporal average of q in the fully-developped regime. The case of externally forced flow is not as straightforward, as starting the forcing at t = 0 (either u or n) triggers nonlinear effects which lead to a change in the flow dynamics. Instead of one, we have two relevant mean flows: a) that associated with the initial attractor characterizing the unforced system at t = 0, and b) that characterizing the final attractor of the forced system as t → ∞. During the transient, the mean flow is not defined. If we want to extend the framework of [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF] to externally forced flows, we therefore need to be very careful about choosing the appropriate reference mean flow. Using the resolvent operator about the mean flow, we can only derive input-output models relevant at the initial and final times. In the open-loop case, we simply have

y = G q n + w, (2.21)
where

G q (s) := C, R q (s)B (2.22)
models the transfer function from n to y and w collects the various remaining terms in the form of an 'output disturbance'. This simple relation is represented by the block diagram in figure 3(b). In the closed-loop case, the transfer function model becomes

G K q (s) := C, R K q (s)B (2.23)
and depends on both q and the controller K. This pseudo-transfer function can also be computed from the expression

G K q = G q 1 -G q K , (2.24)
which corresponds to the linear plant G q in feedback loop with the controller K, as can be seen in figure 3(d). The reference field q either characterizes the mean flow of the initial attractor or the final attractor. The discrete evolution of the reference mean flow, and therefore of the input-output model of the system as we add control action, may be used to design an iterative control strategy as explained below.

Iterative control strategy

In this section, we introduce an iterative control strategy to reach the base flow from the permanent regime of oscillation of an uncontrolled resonating flow. Assume we know how to design a first controller K 1 able to decrease the amplitude of oscillation of the unforced flow, but unable to completely stabilize it, can we design a controller correction K ′ 2 such that the total controller K 2 := K 1 + K ′ 2 is able to damp the oscillations further? In general, from a sequence of m successive controller corrections

K m := 0, if m = 0, K ′ 1 + K ′ 2 + • • • + K ′ m , if m > 0,
(2.25) leading to a permanent oscillatory regime, how can we design the next controller correction K ′ m+1 such that K m+1 = K m + K ′ m+1 leads to even smaller oscillation amplitude? The key idea is to consider an input-output model based on the mean flow q m obtained in the presence of the controller K m . Indeed, as shown in the previous section, the response of the coupled system to a perturbation n of the forcing signal is given by the resolvent operator R Km q m , involving both q m and K m . The corresponding transfer function from n to y is then given by G Km q m . To simplify notations, we denote from now on

G j i := G Kj q i .
(2.26)

At each step m + 1, the plant to control is therefore denoted G m m . The goal is to design a controller correction K ′ m+1 such that the target plant

G m+1 m = G m m 1 -G m m K ′ m+1 (2.27)
has some desired properties, leading to smaller oscillations of the output y. When the additional controller K ′ m+1 is added in feedback loop, the flow transiently evolves until a new dynamical equilibrium is reached, characterized by a new mean state q m+1 . Since the effective controller is now

K m+1 = K m + K ′ m+1 , the transfer function of the new plant is G m+1 m+1 .
The iteration then proceeds until convergence, i.e. until some fixed point is hopefully reached.

The momentum equation with feedback reads

∂ t u + (u • ∇)u = -∇p + Re -1 ∇ 2 u + B u t 0 κ m (τ )y(t -τ ) dτ, (2.28)
where B u is the component of B forcing the momentum equation and κ m is the impulse response of the controller

K m . If a fixed-point q ∞ := (u ∞ , p ∞ ) is reached in the limit t → ∞, then it must satisfy (u ∞ • ∇)u ∞ = -∇p ∞ + Re -1 ∇ 2 u ∞ + B u K m (0)( C, q ∞ -Q + Y ).
(2.29)

In order to reach the base flow, i.e. q ∞ = q b , we must ensure that q b is a solution of (2.29), and that it is the only one. This is the case if the controller has zero static gain K m (0) = 0, and we therefore choose to impose at every step j. In the special case where the sensor measures a deviation with respect to the base flow by satisfying

K ′ j (0) = 0 (2.30) n n n n y + + y + + y + + y + + u (a) (b) (c) (d) w w K ′ 1 u ′ 1 u ′ 1 w w K 2 G 0 0 G 1 1 G 0 1 G 1 1 G 0 1 G 0 2 G 2 2 + + t = t + 1 t = t - 2 t = t + 2 t = t - 3 K ′ 1 K ′ 1 u ′ 2 u ′ 2 K ′ 1 = K 1 K ′ 2 K ′ 2 u u = u ′ 1 u = u ′ 1 G 1 0 G 2 1
Y = C, Q -q b , (2.31)
then q b remains a fixed-point of (2.29) even if K m (0) = 0, but there may be alternative fixed points q ∞ as well. It is therefore crucial to 'disconnect' the controller in the steady regime by imposing (2.30); then q b is for sure the only fixed point for any choice of sensor (C, Y ).

Although the iterative method bears resemblance with the gain scheduling approach of [START_REF] Khalil | Nonlinear systems[END_REF] and [START_REF] Högberg | Relaminarization of Reτ = 100 turbulence using gain scheduling and linear state-feedback control[END_REF], there are fundamental differences which are emphasized later in the discussion, see section §6.2. We stress again that the concept of a transfer function from n to y is ill-defined in the context of nonlinear flows. Our method relies on the assumption that the system responds to external forcing in a LTI fashion when it is close to a state of dynamical equilibrium. The quantity that we will call 'transfer function', and that we will denote G in general, is therefore only defined in two limit cases: either a) the system is in a dynamical equilibrium characterized by q m and K m , or b) the system is out of equilibrium because a controller correction K ′ m+1 has just been added, but the flow has not yet evolved so the mean flow q m remains relevant over some period of time. Denoting t m+1 the instant when the controller correction

K ′ m+1 is added, we have G = G m m at t = t - m+1 and G = G m+1 m at t = t + m+1 . And then we reach G = G m+1 m+1 at t = t - m+2 .
But between the two instants t m+1 and t m+2 , the transfer function model G cannot be defined as there is no mean state relevant to the transiently evolving flow.

Finally, according to equations (2.24) and (2.25), each transfer function G j i can be seen as the interconnection

G j i = G 0 i 1 -G 0 i (K ′ 1 + • • • + K ′ j )
(2.32) between the plant G 0 i characterizing the mean flow alone (even though it cannot be sustained without controllers for i = 0), and the controller corrections K ′ 1 to K ′ j plugged in parallel. The block diagrams at times t = t + 1 , t ± 2 and t - 3 are shown in figure 4.

2.4. Relation between poles of the transfer operators G j i and instability modes In the present context, modal analysis of the Jacobian operator A j i is used to find the pole singularities of the resolvent operator R j i (with obvious notations), and therefore of the transfer function G j i . We are concerned with the solutions of the generalized eigenvalue problem (s l B -A j i (s l ))q ′ l = 0, (2.33) with complex eigenvalues s l = iω l + σ l of growth rate σ l , frequency ω l and eigenmodes q ′ l (the method used to obtain approximate solutions of (2.33) will be explained in §4.1). In our case, resonant frequencies will be associated with specific eigenmodes of the mean-flow q i , coupled with controller K j . These modes cannot be directly interpreted as instability modes if the mean flow is different from the base flow. However, the Jacobian operator A j i may be seen as a perturbation of operator A 0 i with a correction term proportional to K j according to equation (2.16), i.e.

A j i = A 0 i + K j B C, • .
(2.34)

A continuous connection can therefore be established between the eigenvalues of the two operators. The Jacobian operator A 0 i may in turn be seen as a perturbation of the Jacobian operator about the base flow A 0 q b , with a correction term due to mean flow distortion ∆q i , where

∆q := q -q b .
(2.35) Therefore, all instability modes of the base flow may be unambiguously identified with some poles of G j i , thereby providing the latter with a physical meaning. This key advantage of our resolvent-based model will be illustrated in §3.5.

Unforced flow

Numerical model and methods

The control strategy is implemented on a well-known open-cavity flow configuration, which was considered in a series of past papers [START_REF] Barbagallo | Closed-loop control of an open cavity flow using reduced-order models[END_REF][START_REF] Barbagallo | Input-output measures for model reduction and closed-loop control: application to global modes[END_REF][START_REF] Dergham | Model reduction for fluids using frequential snapshots[END_REF][START_REF] Poussot-Vassal | Parametric reduced order dynamical model construction of a fluid flow control problem[END_REF]Schmid & Sipp 2016;Sipp & Schmid 2016) following [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF]. The reader is referred to this initial paper for a comprehensive description of the numerical setup. The cavity of length L = 1 and depth D has an aspect ratio L/D = 1, as shown in figure 1. In the cartesian coordinate system (x 1 , x 2 ), the velocity components are denoted u = (u 1 , u 2 ). With the upstream edge of the cavity defining the origin (0, 0) and the incoming flow aligned with the x 1 -direction, the upstream boundary layer starts forming at x 1 = -0.4 by switching the lower-wall boundary condition from stress-free ∂ 2 u 1 = 0 to no-slip u 1 = 0 (in both cases, the walls are impermeable, i.e. u 2 = 0). The inlet velocity (U = 1, 0) at x 1 = -1.2 is uniform and a standard outflow condition is prescribed at the outlet, which is located at x 1 = 2.5. A stress-free boundary condition is implemented on the upper wall at x 2 = 0.5 and at the downstream-end of the lower wall, from x 1 = 1.75 to x 1 = 2.5. The two-dimensional incompressible Navier-Stokes equations are solved using the freely available finite-element software FreeFem++ (www.freefem.org). The primitive variables (u 1 , u 2 , p) are discretized on a mesh of ∼ 200, 000 (P 1b, P 1b, P 1) Taylor-Hood elements, yielding a total of ∼ 700, 000 degrees of freedom. The time-stepping code, which was made available online (https://github.com/denissipp/AMR Sipp Schmid 2016) as complementary material to the review article by Sipp & Schmid (2016), is semi-implicit, with the nonlinear terms extrapolated with a second order Adams-Bashforth scheme.

We consider the nonlinear dynamics of the linearly unstable flow at a Reynolds number of Re = 7500. The dynamical equilibrium that we consider for our control problem is obtained by initializing the simulation with the base flow q b -which was computed using Newton's iteration (Sipp & Schmid 2016) -plus a small contribution from the most unstable linear mode. After a long transient, characterized in §3.3, the flow reaches a quasiperiodic attractor, characterized in §3.4.

Modal analyses were carried out using the ARPACK implementation of Arnoldi's method [START_REF] Lehoucq | ARPACK users' guide: Solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF], coupled with the multifrontal sparse LU solver MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF]. The shift-invert mode was used to find eigenvalues of the Jacobian operator in the vicinity of a given shift.

Actuator and flow measurements

The upstream-edge actuator is modelled by a Gaussian-shaped volume force term

B = 0, η exp - (x 1 -x 0 1 ) 2 + (x 2 -x 0 2 ) 2 2σ 2 , 0 T , (3.1)
centered around (x 0 1 , x 0 2 ) = (-0.1, 0.02), with a scalar η such that B, B = 1. With the choice of spatial variance σ = 0.0849, the volume force reaches 50% of its maximum value at a radial distance of ≈ 0.1 from (x 0 1 , x 0 2 ), and 1% at ≈ 0.26. In order to characterize the strength of the nonlinearity, we define different measures of the distance between the instantaneous field q and the base flow q b , based on the perturbation field ∆q. The downstream-edge sensor is modelled by a pair (C, Y ) such that

y = 1.1 x1=1 ∂ 2 ∆u 1 | x2=0 dx 1 (3.2)
measures the perturbation wall friction at the downstream edge of the cavity. Note that this sensor satisfies (2.31), hence the base flow is a fixed-point of the closed-loop system, even if the controller has non-zero static gain. However, as discussed in §2.3, we will still need to impose K ′ j (0) = 0 at every iteration, in order to guarantee that there are no alternative fixed points.

In order to build phase portraits, we also consider three x 2 -perturbation velocity probes located within the shear layer at x 2 = 0 and x 1 = 1/4, 1/2, 3/4. Finally, to characterize the system globally, we consider the perturbation kinetic energy

E := ∆q, B∆q q b , Bq b . ( 3

Transient dynamics

Figure 5 shows timeseries of the global and local measures E and y of the uncontrolled flow. We recall that the simulation is initialized with the base flow + a small perturbation parallel to the most unstable base flow eigenmode. There is a short initial phase of linear growth associated with that mode (see inset (a)), before saturating effects of nonlinearity become active. This phase lasts a few convective time units, after which the flow undergoes a long transient evolution of a few hundred time units before a dynamical equilibrium is eventually reached, that we call state 0. This state appears to contain multiple frequencies, as can be seen in insets (b) and (e). Mean flow distortion is caused by nonlinear interactions between these finite-amplitude oscillatory modes. The deviation is significant since the mean value of y is comparable with its oscillation amplitude, confirming strong nonlinearity. The mean perturbation energy E is of the order of 2% (see figure 5(c)), which corresponds to perturbation velocities of the order of 0.1.

Attractor

The resulting mean flow is shown in figure 6(a), together with the base flow in (b), for comparison. We display streamfunction contours on top of filled contours of the total velocity norm. The growth of a boundary layer at the lower wall is apparent from x 1 = -0.4. This boundary layer becomes a shear layer on top of the cavity, which thickens as it is convected towards the downstream edge. The flow recirculates at low speeds within the cavity, as the incoming flow impacts the downstream wall of the cavity. We note two main differences between mean and base flows: (a) the spreading of the shear layer is more pronounced for the mean flow, (b) from inspection of the streamlines, the recirculating flow rate is also stronger in that case. The enhanced flow rate is caused by stronger entrainment due to the unsteady shear layer.

Figure 6(c) shows a vorticity snapshot of state 0 yielding maximal instantaneous shearstress y, compared to the steady vorticity field of the base flow in panel (d). We notice the spatial development of the Kelvin-Helmholtz instability on top of the shear layer, which is responsible for its enhanced spreading in figure 6(a). The maximum shear stress is caused by the impact of a vortical structure at the downstream corner of the cavity. Another vorticity sheet of opposite sign is visible within the cavity, and marks the border of the recirculating zone. The vorticity levels within that sheet are higher in state 0 than in the base flow, which is consistent with the stronger recirculation in the former case. Finally, the average distance between the mean flow and any instantaneous snapshot measured as q ′ , q ′ / q 0 , q 0 ≈ 1.1 × 10 -3 is very low, so the mean flow can be considered as a good approximation of the full flow at any instant.

The dynamical state may be further characterized in the frequency domain by considering the spectrum of Koopman modes. For fully-developped flows, these can be computed using harmonic averages q′ (ω) := lim

T →∞ 1 T T 0 q ′ (t)e -iωt dt (3.4) (Arbabi & Mezić 2017
). The limit is non-zero for Koopman eigenfrequencies only, and yields the associated Koopman modes of the observable q ′ . In practice, the spectrum of amplitude of the Koopman modes may be approximated using DFT on a finite timeseries of perturbation fields. The discrete set of peaks in the DFT spectrum of y correspond to Koopman eigenfrequencies, visible in figure 7(a). There is no continuous spectrum in the signal and the non-zero values between the peaks correspond to noise arising from the estimation by DFT. The precise nature of the attractor is made clear by embedding the dynamics in a 3-dimensional phase space, using the three x 2 -velocity probes in the shear layer. The resulting attractor has a toroidal topology, as can be seen in figure 7(b). Slicing the attractor when trajectories cross the plane ∆u 2 (1/2, 0) = 0 from below yields a Poincaré map in the [∆u 2 (1/4, 0), ∆u 2 (3/4, 0)]-plane, shown in figure 7(c). The resulting closed curve confirms the nature of the attractor as a 2-torus: a quasiperiodic regime with 2 incommensurate 'fundamental' frequencies that we denote ω l and ω h for low-frequency and high-frequency respectively. The discrete set of peaks in the spectrum is generated by all the nonlinear interactions between these two modes:

ω l ω h 2ω h 3ω h 4ω h ω h -ω l ω h + ω l ω h -2ω l 2ω l (b) -0.05 0
{ω k } := {k l ω l + k h ω h | k := (k l , k h ) ∈ Z 2 } (3.5)
Note that the choice of the 'fundamental' frequencies is somewhat arbitrary, as any pair of incommensurate frequencies is appropriate to define the quasiperiodic frequency set (3.5). However, the choice can be physically motivated by selecting two frequencies associated with natural time scales occurring in the flow, as explained below.

Recall that the base flow was initially perturbed with a 'drop' of the leading eigenmode, i.e. that of largest temporal growth rate. This shear-layer mode initially has a natural frequency of ω = 10.9 [START_REF] Barbagallo | Closed-loop control of an open cavity flow using reduced-order models[END_REF]), as can be verified from the spectrogram in figure 8 (the part of the figure for t > t 1 will be discussed later in §5.2). As the instability develops, the frequency shifts to higher values, reaching a value of 11.74 when the flow is fully developed. We use this value to define a high-frequency scale ω h associated with shear-layer oscillations. During the transient, we also observe the development of a lowfrequency peak at ω l ≈ 2.89. The flow structures associated with ω l and ω h are shown in figure 9(a,b). Panel (b) confirms that the Koopman mode at ω h is associated with a shear-layer mode. On the other hand, the low-frequency mode in panel (a) is localized within the cavity. The frequency is low because the recirculating mean flow within the cavity is slow. The two 'fundamental' frequencies ω h and ω l give rise to harmonics and interaction peaks visible in the spectrogram and the spectrum of state 0 in figure 7(a).

Modal analysis

We illustrate in figure 10 the key property of our resolvent-based model explained in §2.4, i.e. that eigenvalues of the mean flow q 0 correspond to poles of the transfer function G 0 0 , which can be related to instability modes of the base flow q b . On the top panels, we plot the gain of the frequency response associated with the unforced mean flow G 0 0 (iω) (panel (a)) and base flow G b (iω) (panel (b)), the latter being defined as

G b := C, R q b B .
(3.6) c,d) calculated by modal analysis. Eigenvalues are coloured according to criterion (3.7) quantifying the relative contribution of each pole to the input-output relation. The (almost) one-to-one relation between resonance peaks and specific poles is highlighted with dashed vertical lines. These resonant poles can be grouped into a single family for the base flow, highlighted with a grey line in panel (d). This family contains the four unstable modes of the base flow, labelled with greek characters α, β, γ and δ in decreasing order of growth rate and encircled in black. This family of resonant modes can also be identified, by continuity, in the mean flow spectrum. A distinct family of modes a-f , also highlighted with a grey line in panel (c), is responsible for the low-frequency peaks present in |G 0 0 |.

In both cases, we notice the presence of peaks indicating strong receptivity to external forcing on a discrete set of frequencies. In the case of panel (a), these peaks correspond to the frequencies of the Koopman modes. Hence, the unforced flow is most receptive to external forcing at its intrinsic oscillation frequencies, which is an indication that the mean flow model captures relevant properties of the nonlinear system. The peaks of G b , on the other hand, do not match with nonlinear frequencies. In the two cases though, the peaks are associated with resonant poles of the respective model, corresponding to meanflow or base-flow eigenvalues. We therefore plot the spectra of the associated Jacobian operators in panels (c) and (d), below the frequency responses. In order to identify which eigenvalues are responsible for the receptivity peaks, we color the spectra by the criterion

Γ j = | C, q ′ j || p ′ j , B | |σ j | , (3.7)
quantifying the contribution of each mode to the input-output relation [START_REF] Antoulas | Approximation of Large-Scale Dynamical Systems[END_REF][START_REF] Bagheri | Input-output analysis, model reduction and control of the flat-plate boundary layer[END_REF][START_REF] Barbagallo | Closed-loop control of an open cavity flow using reduced-order models[END_REF]. The criterion takes into account the distance 1/|σ l | of each eigenvalue to the imaginary axis, the observability of the direct mode | C, q ′ l | and the controllability of the adjoint mode | p ′ l , B |. We recall that the adjoint modes (p ′ l ) are solutions of the eigenvalue problem (µ l B -A j i (µ l ))p ′ l = 0, where A j i is the adjoint of A j i with respect to the inner product (2.10). Eigenvalues of the adjoint operator are complex conjugate of the eigenvalues (s l ) of the direct operator, i.e. µ l = s * l , and adjoint modes satisfy the biorthogonality relation p ′ k , q ′ l = δ kl with the normalized direct modes (q ′ l ). Using the criterion (3.7), we find two families of modes responsible for the resonance peaks of G 0 0 , and a single family of modes responsible for the peaks of G b . The modes are highlighted with grey lines in panel (c) and (d), intended to guide the eye. Four poles, labelled α, β, γ, δ are responsible for the peaks at ω > 5 in G 0 0 , corresponding to oscillations at ω h , ω h ± ω l and ω h + 2ω l . The other family of mean-flow modes a, b, c, d, e, f are responsible for the low frequency peaks at ω < 5, among which the peak at ω l , coinciding with mode b. The five resonance peaks of the base flow are generated by a single family of 5 poles. From the time-frequency analysis, we know that the high frequency ω h is connected to the most unstable mode of the base flow. Therefore, mode α is the mean-flow analog of that specific base flow eigenmode. The other meanflow modes β, γ, δ may then be identified with the remaining unstable modes of the base flow, hence the labelling in figure 10(c,d).

Conversely, the mean flow mode b responsible for the receptivity peak at ω l in G 0 0 does not have an obvious base flow analog. In fact, the emergence of mode b in the meanflow spectrum seems to arise from a purely nonlinear mechanism comparable to 'lock-in'. Indeed, we notice that the frequency difference between each pair of modes within the set {α, β, γ, δ} is a multiple of the frequency ω l associated with mode b. Distortion of the base flow by the nonlinearities allows these four modes to synchronize with the low-frequency mode b, resulting in a quasiperiodic attractor. A similar observation has already been made regarding low-frequency oscillations of large laminar separation bubbles behind a two-dimensional bump [START_REF] Cherubini | The effects of non-normality and nonlinearity of the Navier-Stokes operator on the dynamics of a large laminar separation bubble[END_REF].

In figure 9, we compare the structure of the mean-flow eigenmodes α and b (panels (c, d)) with the Koopman modes of the unforced flow at ω h and ω l (panels (a, b)). We notice striking similarities between eigenmodes and Koopman modes, in line with observations from past studies recalled in the introduction. In the literature though, the success of mean flow stability analysis is often associated with the presence of marginal eigenvalues [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF][START_REF] Mantič-Lugo | Self-consistent model for the saturation mechanism of the response to harmonic forcing in the backward-facing step flow[END_REF], the so-called RZIF property (real-zero imaginary-frequency) in [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF]. We highlight the fact that the RZIF property is not verified here: modes b, α, β, γ and δ have non-zero, positive growth rates, yet they are associated with the correct frequencies and Koopman modes. This was already noticed by Mettot et al. (2014a) and [START_REF] Meliga | Harmonics generation and the mechanics of saturation in flow over an open cavity: a second-order self-consistent description[END_REF] in opencavity flows (although in the former paper the mean flow was obtained as a fixed-point of a RANS model, not a temporal average of an unsteady RANS simulation). Here, the approximation works because the mean flow yields an excellent approximation of the full flow at any instant [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF], as previously said in §3.4.

Finally, we observe in figure 10 receptivity peaks at low frequencies different from ω l , due to the presence of poles c, a and f close to the imaginary axis. The corresponding eigenmodes are localized within the cavity and resemble mode b. Despite strong receptivity to external forcing, there is no signature of the associated frequencies in the spectrum of the unforced flow. To understand this, we need to consider the endogenous forcing term

f ′ := Ψ q q ′ (3.8)
due to the nonlinear feedback operator. Since the operator Ψ q is time-invariant, its output f ′ is quasi-periodic with the same discrete set of frequencies as its input q ′ , i.e. {ω k } defined in (3.5). More precisely, we can introduce the sets of Koopman modes

{q ′ k } := {q ′ (ω k )}, { f ′ k } := { f ′ (ω k )} (3.9)
respectively associated with q ′ and f ′ . In the fully-developped regime, the governing equations (2.2) can be recast in the form of a harmonic balance relation [START_REF] Khalil | Nonlinear systems[END_REF]) (3.10) between these two sets of Koopman modes (harmonic averages are zero otherwise). Clearly, strong receptivity at frequencies outside the set {ω k } has no impact on the intrinsic dynamics because the endogenous forcings are zero at these frequencies. In short, the resolvent operator displays strong receptivity to external forcing at the natural oscillation frequencies of the unforced flow, but not exclusively.

q′ k = R q (iω k ) f ′ k ,

Designing controller K ′ 1

We explain here the methodology used to synthesize controller K ′

1 , but all subsequent controller corrections are designed using the same methodology. There are two consecutive steps: first model reduction from frequency domain data G 0 0 (iω), then structured H ∞ -synthesis using the ROM.

Model reduction

We convert the frequency response G 0 0 in the form of a state-space model ( A, B, C, D). Although many points O(100) are typically needed to capture the fine details of the frequency response, we seek a state-space model approximating G 0 0 with a lower order O(10) to focus on the essential characteristics of the system that we would like to manipulate. A strictly proper (i.e. D = 0) reduced-order model (ROM) of order 16 is obtained using subspace identification methods [START_REF] Liu | Frequency domain structural system identification by observability range space extraction[END_REF][START_REF] Mckelvey | Subspace-based multivariable system identification from frequency response data[END_REF]. For subsequent iterations, the order of the ROM varies between 16 and 32. The frequency response of the ROM is given by

G 0 0 (iω) = C(iωI -A) -1 B (4.1)
and the associated poles correspond to eigenvalues of A. In general, all quantities related to ROMs will be denoted with a tilde. In figure 11(a,b), we check that the frequency response G 0 0 (iω) (black dotted curve) coincides with the original model G 0 0 (solid red curve) over the frequency range 5 < ω < 20. This range contains the main resonance peaks at ω h , ω h ± ω l and ω h + 2ω l , but not the low frequency ones at ω < 5 which we choose to ignore in order to ensure sufficient order reduction. In panel (c), we check that the poles α, β, γ and δ corresponding to the main resonance peaks are correctly approximated by some eigenvalues of A.

We note in passing the monotonically decreasing phase (except for 'upward jumps' corresponding to unstable modes) with a nearly constant negative slope in panel (b), which characterizes a time-delay between the actuator and the sensor. This time-delay corresponds to advection of the perturbations over the cavity of length L = 1 at the mean velocity κU in the shear layer. Fitting the data with a straight line in the range range 5 ω 25 leads to a value of κ ≈ 0.73, consistent with the mean flow in figure 6(a). This value of κ is larger than the value 0.53 obtained by [START_REF] Barbagallo | Closed-loop control of an open cavity flow using reduced-order models[END_REF] for the advection of pressure perturbations by the base flow at x 2 = 0 above the cavity. This difference is mainly due to the fact that the forcing field B is localized within the incoming boundary layer at x 0 2 = 0.02 > 0 and not directly at the wall.

4.2. Structured H ∞ -synthesis After designing a ROM, we seek a controller K ′ 1 outputting the forcing signal

u ′ 1 := K ′ 1 y (4.2)
from the measurement y. When the plant G 0 0 is interconnected with K ′ 1 , the closed-loop is characterized by four transfer functions

y u ′ 1 = G 0 0 S S T K ′ 1 S n w , (4.3)
linking the two inputs n and w to the two outputs y and u ′ 1 . These four transfers are combinations of only two interconnected blocks G 0 0 and K ′ 1 , therefore they are not independent of each other. In particular, the sensitivity function S and complementary sensitivity function T , respectively defined as

S = 1 1 -G 0 0 K ′ 1 and T = G 0 0 K ′ 1 1 -G 0 0 K ′ 1 , (4.4)
satisfy the relation S -T = 1. The controller K ′ 1 is designed in order to achieve desired specifications on the four closed-loop transfer functions. But because the four transfer functions are not independent of each other, we can only constrain each of them in a limited frequency range. Four filters W S , W T , W GS , W K ′ S are therefore designed in order to shape each closed-loop transfer function adequately and achieve specific design goals to be listed in the next subsections. The procedure amounts to minimizing a positive constant ρ bounding the H ∞ -norm of the four weighted transfer functions

W S S ∞ ρ, (4.5) W T T ∞ ρ, (4.6) W GS G 0 0 S ∞ ρ, (4.7) W K ′ S K ′ 1 S ∞ ρ, (4.8)
while simultaneously ensuring that the closed-loop transfers and the controller K ′ 1 remain stable. We recall that for a single-input single-output (SISO) system, the H ∞ -norm of a stable transfer function H is given by the maximum of the modulus of the frequency response for all real values of ω. The design of controller K ′ 1 therefore rests on the definition of four inequality constraints on the closed-loop transfer functions. The constraints are enforced by designing well-chosen frequency domain templates 1/W S , 1/W T , 1/W GS and 1/W KS . These are defined as filters depending on free parameters which are calibrated manually before ρ is minimized. If one of the closed-loop constraints is not met, the free parameters are adjusted and ρ optimized again, in an iterative fashion until all constraints are satisfied.

H ∞ = max
The present approach targets the shaping of all four closed-loop transfer functions and differs from the mixed-sensitivity method, which only targets S, T and sometimes K ′ 1 S (Henning & King 2007;Henning et al. 2007;[START_REF] Williams | Unsteady lift suppression with a robust closed loop controller[END_REF], but not G 0 0 S. The advantage of the four-block method over mixed-sensitivity will become clear in §4.2.1. We note that directly shaping the closed-loop transfer functions is a much more challenging task than shaping the loop gain G 0 0 K ′ 1 , as done in the loop-shaping approach [START_REF] Dahan | Feedback control for form-drag reduction on a bluff body with a blunt trailing edge[END_REF][START_REF] Jones | Modelling for robust feedback control of fluid flows[END_REF][START_REF] Dalla Longa | Reducing the pressure drag of a d-shaped bluff body using linear feedback control[END_REF][START_REF] Li | Feedback control of combustion instabilities from within limit cycle oscillations using H∞-infinity loop-shaping and the ν-gap metric[END_REF]. Loop shaping takes advantage of the linearity of the loop gain with respect to the controller K ′ 1 in order to indirectly achieve constraints on the closed-loop at low cost, whereas both the mixed-sensitivity and the four-block frameworks target the closed-loop in a more direct fashion.

Here we use a structured H ∞ framework, which means that the controller order may be chosen a priori, i.e. independently of the order of the ROM. In our case, we fix the order of K ′ 1 and all subsequent controller corrections to 10. This is another key difference with the loop-shaping approach, where the order of the controller cannot be chosen freely and is always larger than that of the ROM. However, computing fixed-order controllers for the H ∞ -synthesis problem is known to be NP-hard [START_REF] Blondel | NP-hardness of some linear control design problems[END_REF] and one has to rely on local optimization techniques to compute simple and practical controllers. In this work, we have used the control software hinfstruct from the MATLAB Control Toolbox, which is based on a non-convex non-smooth bundle technique developed in [START_REF] Apkarian | Nonsmooth H∞ synthesis[END_REF].

In the rest of this section, we will explain the specifications on the closed-loop transfer functions and the choice of appropriate frequency templates. The final transfer functions (black lines) and frequency templates (blue lines) are shown in figure 12, before (solid) and after (dashed) optimization of the parameter ρ to a final value of 0.6. 4.2.1. Shaping G 0 0 S: attenuation of the main resonance peak The goal of the controller is to suppress intrinsic resonances. To 'desensitize' the flow to input perturbations n at the main resonance frequency ω = ω h of G 0 0 (iω), we must ensure that

|S(iω h )| < 1, (4.14) such that the closed-loop gain | G 1 0 | = | G 0 0 S| (4.15)
be smaller than the open-loop gain | G 0 0 | at this frequency. This is achieved by designing a 'bell-shaped' frequency template

1/W GS = 1/k GS s 2 /α 2 + 2ξω h s + ω 2 h /α 2 s 2 + 2ξω h s + ω 2 h (4.16)
with α 2 > 1 (see figure 12(a)) applied to the closed-loop transfer function G 0 0 S. The gain |1/W GS (iω)| varies between a minimum value of 1/(α 2 k GS ), reached asymptotically when ω → 0 and ω → ∞, and a maximum value of 1/k GS , reached at the resonance frequency ω = ω h . The parameter ξ determines the width of the filter. By choosing Since the optimal value of ρ is 0.6 in our case, we enforce a reduction of 60% in the height of the main resonance peak at ω h . Shaping G 0 0 S in order to constrain S may seem counterintuitive at first sight: why not directly shaping S? The aim of this 'trick' is to prevent pole cancellation by a zero of the controller. Indeed, cancelling the pole α with a zero of K ′ 1 is not a robust way to suppress the resonance at ω h since the mean flow model is only an approximation of the true plant. Worse still, the unstable pole α will remain present in the closed-loop through the transfer G 0 0 S. A much better alternative is to strongly damp the growth rate associated with this pole by shaping G 0 0 S instead of S. When targetting the sensitivity function, as in the mixed-sensitivity approach, pole cancellation would typically occur since S is a function of the product G 0 0 K ′ 1 , i.e. the loop gain, only. The same problem occurs in the loopshaping framework, which directly targets the loop gain. When targetting G 0 0 S instead, pole cancellation does not occur since that transfer depends on both the loop gain and G 0 0 . The four-block approach is therefore well-suited for the control of unstable plants.

k GS = 1/| G 0 0 (iω h )| (4.17
|K ′ 1 S| ω K ′ 1 S |T | ω T 1/|W T | ρ/|W T | | G 0 0 S| ω G 0 0 G 0 0 S 1/|W G 0 0 S | ρ/|W G 0 0 S | |S| ω S 1/|W S | ρ/|W S |
In fact, with appropriate choices for α and ξ, we see in figure 12(a) that we manage to suppress all three resonance peaks at ω h and ω h ± ω l , associated with modes α, β and γ.

As a result, we verify in panel (b) that the gain of the sensitivity function is less than 1 at these three frequencies. We also verify in this panel that ρ = |S(iω h )| = 0.6. We will later verify in figure 17 from §5.3 that the corresponding poles α, β and γ are indeed damped instead of being cancelled. 4.2.2. Shaping S: stability robustness (modulus margin)

The poles of all four closed-loop transfer functions are given by the zeros of 1 -G 0 0 K ′ 1 . To ensure closed-loop stability, we must therefore ensure that all zeros lie in the appropriate half-plane. To guarantee stability robustness, we also ask that these zeros be sufficiently far from the imaginary axis, where the closed-loop becomes marginally stable. Equivalently, for a stable closed-loop, we may ask that the modulus margin (not to be confused with the gain margin)

∆M = min ω∈R |1 -G 0 0 (iω)K ′ 1 (iω)| (4.19)
be sufficiently large. In practice, we ask that ∆M > -6dB, or ∆M > 0.5 in linear scale, which amounts to bounding the sensitivity function

S ∞ < 2. (4.20)
The constraint on S(iω) being independent of ω, we consider a constant frequency template 1/W S and optimize ρ such that

|S(iω)| ρ |W S | < 2 (4.21)
for any real ω. In figure 12(b), constraint (4.20) is represented by the shaded region. We see that the value of ρ = 0.6 reached after optimization did not allow the constraint to be satisfied in the vicinity of ω ≈ 13, with the chosen W S . However, the sensitivity gain being almost equal to 2 there, the resulting controller is deemed satisfactory without further changing W S .

Shaping T : enforcing zero controller static gain

To ensure that the base flow is the only fixed-point of the closed-loop, we need to ensure the controller has zero static gain, as discussed in §2.3. Imposing (2.30), i.e. K ′ 1 (0) = 0, is equivalent to a condition on the complementary sensitivity function

T (0) = 0, (4.22)
which can be approximately achieved using constraint (4.11) and a filter W T with a sufficiently large static gain

|W T (0)| ≫ 1, such that |T (0)| < ρ |W T (0)| ≪ 1. (4.23)
The frequency template 1/W T is therefore defined as a high-pass filter, as can be verified in figure 12(c).

4.2.4. Shaping K ′ 1 S: performance robustness to modelling errors at high frequency Finally, the remaining transfer function K ′ 1 S may be shaped in order to 'desensitize' the controller at high frequency. This can be useful if modelling errors are large at high frequencies and if the signal y is stochastic or contaminated with noise. This would typically be the case in an amplifier flow or in the context of turbulence. To achieve sufficient 'roll-off' of the controller gain at high frequency, a low-pass frequency template 1/W K ′ S may be used. In our case, we expect our model to become inaccurate at large frequencies since the resolvent operator about the mean flow only captures spatially coherent structures typically developing at low to moderate frequencies. However, our signal y being deterministic, non-chaotic, noise-free and therefore weakly energetic at high frequencies (cf. spectrum in figure 7(a)), imposing steep roll-off at high-frequency is unnecessary. We therefore choose 1/W K ′ S as a very large constant here, and simply verify a posteriori in figure 12(d) that the obtained controller is strictly proper, i.e. that |K ′ 1 (iω)| → 0 as ω → ∞.

Controlled flow

Timescales

A sequence of five robust controllers was found to lead the flow from state 0 to the base flow. When a controller correction K ′ m is added to the loop, its internal state is initialized to 0. The internal states of the other controllers K ′ j<m already in the loop are not reset to 0, hence the control signal u evolves continuously over time. Figure 13 shows the full time-evolution of the system as the sequence of controllers is applied. In this plot, we concentrate on the transient evolution from one state to the next. The origin of time is reset to zero at each iteration m, by subtracting the start-up time t m of a new controller K ′ m . The time-axis is in log scale, and different shades of grey are used to highlight the different time scales involved in the journey of the system from one state to the next. When a new controller is switched on, the flow rapidly evolves over a short timescale of the order of O(1), which seems sufficient for the system to latch onto a new state. This timescale is indicated with a white background. Subsequent evolution of the local measurement y is insignificant, but the perturbation kinetic energy keeps evolving considerably over the next O(100) time units. This slow relaxation is indicated with a light grey background, and a dark grey background is used to denote convergence to a new equilibrium. The separation of timescales in the transient evolution is once again related to the two physical processes occuring in the flow: fast oscillations of the shear layer and slow recirculation within the cavity. Indeed, due to its position, the local measurement y is mostly sensitive to the former process and very little to the latter. Conversely, the global measurement E is, by definition, sensitive to both. The control signal is characterized by a third, ultra-fast, timescale O(0.1), highlighted in very dark grey and visible at iteration 1 only. This timescale is intrinsic to the controller and the large 'overshoot' observed at the first iteration is due to the transition from a two-torus to a limit-cycle with a different fundamental frequency. Such qualitative transition does not occur in subsequent iterations, as we shall see next.

In the last iteration, all three quantities vanish as the system is attracted towards the base flow, which is a solution of the Navier-Stokes equation without forcing. At this point, the power consumption of the controllers is infinitesimal (or proportional to incoming noise in a realistic case) and the cavity drag is minimal, so the control strategy is profitable in the long run, regardless of the initial cost. The exponential decay of all quantities is governed by the least damped pole of the last transfer function G 5 5 .

Attractors

In figure 8 from §3, we plotted a spectrogram of the timeseries y(t), which shows the evolution of the frequency content of the flow as it evolves. We comment here the part for t > t 1 concerning the controlled flow. We immediately notice in the spectrogram that there is less and less energy at high frequencies as we iterate, which indicates that the nonlinearity becomes weaker and weaker, before a linear regime is eventually reached at iteration 5. Initially, the dynamics is on a two-torus characterized by the two fundamental frequencies ω l and ω h . But soon after t 1 , the flow transitions to a limit cycle with a fundamental frequency close to the previous harmonic ω h + ω l . During subsequent iterations, the flow stays on (nearly) that same attractor, but the fundamental frequency slightly decreases. In fact, in the second iteration (between t 2 and t 3 ), the attractor is not exactly a pure limit cycle but an almost degenerate two-torus, which is strongly dominated by its fast fundamental frequency. The trace of interactions between the highfrequency fundamental and a non-zero low-frequency mode at ω ≈ 2.5 is visible soon after t 2 in figure 8.

Focussing on the permanent states in figure 13, we notice that each controller leads to a reduction of up to a decade in the perturbation kinetic energy (panel (a)). The decrease of the oscillation amplitude at the sensor y is less dramatic (panel (b)). In figure 14, we project each attractor in the phase plane [∆u 2 (1/4, 0), ∆u 2 (3/4, 0)]. We observe a qualitative change in the topology of the attractor from state 0 (two-dimensional) to state 1 (one-dimensional). After that, each additional controller further reduces the amplitude of shear-layer oscillations, without further qualitative change in the dynamics.

In figure 15, we also plot vorticity snapshots for each state, at a time of maximal shear stress y. The reduction of the mean recirculation and the suppression of the low-frequency cavity mode are responsible for the strong weakening of the vorticity sheet within the cavity between states 0 and 1. The snapshots also illustrate the weakening of the spatial instability in the shear-layer: roll-up is delayed and the resulting vortices are therefore weaker as they impact the downstream edge of the cavity. We also notice a reduction in the size of the vortices between states 0 and 1, which is consistent with the increase in the oscillation frequency observed in figure 8. m-1 (dashed lines) at t = t + m to G m m (solid lines) at t = t - m+1 (with t6 → ∞), under the effects of nonlinearity. Same labels as figure 15. Vertical lines indicate Koopman frequencies of the permanent state obtained with the controller Km (no unsteady peak for m = 5, because the system is steady). We also indicate which poles of G m m are responsible for each receptivity peak. Mode labels are coloured when their receptivity peak matches with a Koopman frequency of the flow. controller correction K ′ m+1 is added to the loop, in order to desensitize the plant with respect to incoming perturbations occurring at resonant frequencies. As previously said in §4.2.1, each controller targets the receptivity peaks matching with Koopman modes (indicated with solid lines), using a 'bell-shaped' filter (cf. equation (4.16)) centered about the dominant Koopman mode. This procedure leads to the shaped plant G m+1 m (dashed lines). The flow being out of equilibrium, it transiently evolves until a new dynamical state is reached. The new plant G is then characterized by the transfer function G m+1 m+1 (solid lines) and the iteration proceeds. In the last step, the flow state 4 is very close to the base flow already and the designed closed-loop G 4 5 is therefore nearly equal to the effective plant G 5 5 .

In figure 17, we plot the poles of each transfer function G during the iteration. The poles of G m m are indicated with solid symbols and those of the shaped plant G m+1 m are indicated with open symbols. In fact, the poles are computed from a ROM G obtained by subspace identification of G, at each step. We only follow the evolution of the four poles α, β, γ, δ captured by the ROM. We checked in §4.1 that the poles are well approximated in the case of G 0 0 (see figure 11(c)). The plant G m m always has unstable poles, until the very last iteration. During the 'shaping' phase (dashed arrows) G m m → G m+1 m , the plant is stabilized. The growth rate of all four resonant poles is damped (except pole γ in the last iteration; panel (f )). Pushing the poles far from the imaginary axis attenuates the receptivity peaks in the transfer function, as observed in figure 16 and previously explained in §4.2.1. The most strongly damped pole is associated to the frequency of the dominant Koopman mode, specifically targeted by the 'bell-shaped' template. Then, during the relaxation phase G m+1 m → G m+1 m+1 , all the poles drift upwards again (solid arrows) and settle closer to the imaginary axis, causing the resurgence of sharp peaks in the next plant G m+1 m+1 . In the last step, the flow is weakly nonlinear and the drift of the poles during the relaxation phase is sufficiently small for all of them to remain in the stable half-plane, which leads to convergence to the base flow.

We finally highlight the fact that, at each step m, the fundamental oscillation frequency of the fully-developped flow matches with a specific pole of G m m . In the case G 0 0 of the unforced flow, there is even a correspondence between five Koopman frequencies and five poles. This remarkable observation confirms that each linear model based on the mean state captures the intrinsic dynamics of the corresponding nonlinear system {flow + controller K m }. Given the fact that each pole may be associated with a base flow eigenmode, we can identify which instability modes are driving the observed nonlinear branches. The association is made explicit in the spectrogram of figure 8, where mode labels are indicated on top of each frequency peak. Initially, the flow is in a quasiperiodic regime due to an interaction between all four modes α, β, γ and δ, which synchronize with mode b. Then the controlled flow switches to a limit-cycle driven by mode β, while the signature of the other modes is completely suppressed from the output signal y. In the unforced case, the leading Koopman mode is associated with mode α, whereas it is associated with mode β when K ′ 1 is in the loop, which explains the sudden frequency jump in the spectrogram. Interestingly, we note the complete disappearance of the low frequency at ω l , despite the fact that controller K ′ 1 did not target the corresponding mode b. The suppression of this oscillation probably occurred through a nonlinear mechanism: by damping modes β, γ and δ corresponding to the oscillations at ω h -ω l , ω h + ω l and ω h + 2ω l , the controller K ′ 1 has indirectly targetted the low frequency ω l . In subsequent iterations m > 1, the flow remains on a limit-cycle governed by mode β, which explains the smooth evolution of the dominant frequency in the spectrogram (figure 8), and the qualitatively similar phase portraits in figure 14.

Discussion

Alternative route to the base flow and robustness of the

We checked the robustness of our method by using a different set of controllers for the last 3 iterations. This choice led to a different route in phase space (not shown), exciting modes γ and δ and suppressing the oscillation initially associated with mode β altogether. It was however still possible to reach the base flow in five steps, which shows that different sequences of controllers may be equally suitable.

Intuitively, the number of iterations increases with the nonlinearity and decreases with

m m-1 = G 0 m-1 /(1 -G 0 m-1 K m ) at step m (cf. equation (2.32)).
To converge in m iterations, we therefore need G 0 m-1 sufficiently close to G b , i.e. a weakly nonlinear state m -1, and a sufficiently robust controller K m . This is exactly what happens for us after m = 5 iterations. With stronger initial nonlinearity, a higher number of iterations may be required, but provided each intermediate input-output model is able to correctly capture the intrinsic dynamics of the system, one would still expect convergence to the base flow, although there is no theoretical proof and the method remains heuristic at this point. 6.2. Key differences with a gain scheduling approach By updating the controller at each iteration according to the evolution of the mean flow, the method is reminiscent of the gain scheduling approach [START_REF] Khalil | Nonlinear systems[END_REF][START_REF] Högberg | Relaminarization of Reτ = 100 turbulence using gain scheduling and linear state-feedback control[END_REF]. There are however key differences listed below.

First, we do not know in advance what mean flows will be obtained, so rather than being scheduled in advance, the controller gain is adapted during the iteration, depending on the new mean flow emerging at each step. In [START_REF] Högberg | Relaminarization of Reτ = 100 turbulence using gain scheduling and linear state-feedback control[END_REF], a sequence of mean flows was determined a priori, by interpolating between the fully-developped turbulent state and the laminar state. The mean flow models therefore characterized the transient evolution of the globally stable (channel) flow during relaminarization. In our case, the mean flows characterize different permanent states visited during the iteration. The attractors may be driven by different instability modes, resulting in qualitatively different mean flows which cannot be predicted by interpolation.

Second, instead of designing a new controller K m+1 from the sole knowledge of the mean flow m, we design a controller correction K ′ m+1 = K m+1 -K m based on the knowledge of both the mean flow m and the controller K m . Again, we are not following the transient evolution of a globally stable flow, but a sequence of permanent states visited by an unstable flow coupled with some controllers. The plant model G 0 m based on the mean flow m alone does not caracterize the response to a small change in actuation of the coupled system, and the adequate transfer function is G m m . As can be seen in figure 18 for m = 5, the transfer functions G 0 m and G m m may be quite different indeed. Finally, in [START_REF] Högberg | Relaminarization of Reτ = 100 turbulence using gain scheduling and linear state-feedback control[END_REF], the order of the controller remains fixed during the iteration whereas that of the total controller K m = m j=1 K ′ j increases in our case. This drawback of our approach may not be an issue for implementation, since each control law u ′ j may be computed independently, and then added to form u m = m j=1 u ′ j . Since each K ′ j has the same order, the computational cost increases only linearly with the number of iterations m. The computational time may even be kept constant using parallel computing. Alternatively, a low-order model may be obtained from the total controller K m at each step. In our case, we successfully reduced the number of states of K 5 from 50 to 12 with no loss in performance (using the modred command from Matlab, which is based on balanced truncation).

Conclusion

The present study proposes a novel approach for the control of strongly nonlinear resonating flows using linear methods only. The approach is illustrated on the case of open cavity flow at Re = 7500, i.e. far from the linear instability threshold which is at Re c = 4140. Starting from a quasiperiodic unactuated regime, we manage to completely suppress the oscillations with a sequence of five robust controllers. At the end of the process, the base flow is reached, with stabilized linear dynamics. The method relies on two key ideas:

(i) modelling of the system as an LTI plant using the resolvent operator about the mean flow, (ii) a Newton's-type iterative procedure based on the successive mean flows.

The first point (i) builds on previous observations that linear analysis about the mean state yields quantitative information about the underlying nonlinear dynamical equilibrium. Modal analysis is indeed able to capture the structure and frequencies of self-sustained oscillations, while resolvent analysis successfully captures receptivity to forcing, even in turbulent flows. Using the resolvent operator in order to derive a sequence of input-output models, we show that it takes no more than a series of time-averaged flows to design a successful feedback loop. Not only is the method simple, general and computationally inexpensive, it is also physically insightful. Furthermore, LTI models of nonlinear flows derived from the mean state are also unambiguously defined, as opposed to linear models identified from input-output data. This is particularly useful in the case of unstable base flows, which inevitably lead to nonlinear dynamics even in the most quiet environment. When the base flow is reached, the system is truly linear and the proposed definition matches exactly with the actual definition of the transfer function. In essence, the framework initially introduced by McKeon & Sharma (2010) naturally lends itself to the modelling of an LTI input-output relation for nonlinear flows. Although seemingly optimal in some sense, there is however no proof that the resolvent operator about the mean flow is the best option, or work in every cases, and these questions will be addressed in future work.

The second point (ii) also appears to be fairly novel, although it is in principle independent of (i), which means that an iterative procedure could also be carried out with a different modelling strategy for defining an LTI model at each step. The key idea here is to resort to a sequence of linear approximations in order to solve a nonlinear control problem. For sufficiently strong nonlinearity, designing a single low-dimensional LTI controller able to destroy the initial attractor and drive the system directly to the fixed point is a challenge. Instead of designing complex nonlinear models, or crafting elaborate nonlinear control laws, we propose to use a simple sequence of linear controllers. Each controller disrupts the current dynamical equilibrium and drives the system to a new one, closer to the fixed point in phase space. After a few iterations, the coupled system becomes attracted to the fixed point and flow oscillations are fully suppressed. Note that although only linear controllers are used, the entire procedure is nonlinear in essence, as the effect of each controller depends fundamentally on the choice of the past controllers. And although the final controller stabilizes the base flow once reached, it is not guaranteed to drive the system from an initial condition on the attractor to the base flow. The method set forth in this paper is in some sense analog to Newton's iteration: a series of linear models and controllers are derived, until fluctuations arising from nonlinearity are fully cancelled.

From a practical perspective, the proposed method has one clear advantage: it only requires knowledge of the mean flow, at every iteration. There is no need for an adjoint simulation either, and, in the case of turbulent flows, the method might be applicable from a time-averaged unsteady RANS simulation. Computing the transfer function requires the formation of the Jacobian matrix, followed by the resolution of a linear system for each frequency. For two-dimensional, or equivalently three-dimensional problems with one invariant direction, this is computationally inexpensive. For fully threedimensional mean flows, matrix-free methods may be considered instead. A disadvantage of the method though, is that it requires knowledge of the full mean state. However, data assimilation is an active field of research in fluid mechanics, and techniques have recently been developped for two-dimensional mean-flow reconstruction from underresolved/incomplete/noisy measurements [START_REF] Foures | A data-assimilation method for reynolds-averaged navier-stokes-driven mean flow reconstruction[END_REF][START_REF] Symon | Data assimilation of mean velocity from 2d piv measurements of flow over an idealized airfoil[END_REF]. For the method to be applicable outside the lab though, methods should be developped to infer the mean flow from localized sensors only.

Future work will be dedicated to understanding the rationale behind the success of mean-flow based input-output models and their limitations. It would also be interesting to analyze the convergence properties of the iterative method from a theoretical viewpoint. Applicability of the method to noise amplifiers, chaotic flows and turbulent flows also remains to be investigated.

Figure 1 .

 1 Figure 1. Numerical model of open cavity flow, with upstream (volume force) actuator B (in yellow) and downstream (shear-stress) sensor (C, Y ) (in blue). The shear-layer instability and cavity recirculation are symbolized with red lines. There are three x2-velocity probes in the shear layer at x2 = 0 and respectively x1 = 1/4, 1/2, 3/4. Grey thick lines indicate no-slip boundary conditions.

Figure 3 .

 3 Figure 3. Block diagram of the forced Navier-Stokes equations in the (a) open-loop case and the (c) closed-loop case. Equivalent representation with an open-loop transfer function GQ from n to y and an output disturbance w in (b) the open-loop and (d) the closed-loop.

Figure 4 .

 4 Figure 4. Evolution of the block diagram representation of the system during the iterative procedure: (a) t = t + 1 , K ′ 1 has juste been added so the flow alone is still characterized by the plant G 0 0 ; (b) t = t - 2 the feedback loop between the fluid and the controller K1 := K ′ 1 has converged to a new mean state; (c) t = t + 2 , K ′ 2 has juste been added so the flow alone is still characterized by the plant G 0 1 ; (d) t = t - 2 the feedback loop between the fluid and the controller K2 := K ′ 1 + K ′ 2 has converged to a new mean state.

Figure 5 .

 5 Figure 5. Timeseries of (a, b, c) perturbation kinetic energy E and (d, e, f ) shear-stress measurement y of the unforced flow, with sliding mean shown as a black dotted line. Zooms on the initial and final stages are shown in insets.

Figure 6 .

 6 Figure 6. Comparison between (a, c) state 0 and (b, d) the base flow. (a, b) Velocity 2-norm of (a) uncontrolled mean flow u0 (computed by averaging over 200 time units, after the transient) and (b) base flow u b , with associated contours of the streamfunction ψ. The streamfunction is defined by the relations ∂2ψ = v1, ∂1ψ = -v2 and ψ = 0 on the lower wall (solid contours for positive increments δψ = 0.05 from ψ = 0; dashed contours for negative increments δψ = -0.01 from ψ = 0). (c, d) Vorticity fields: (c) snapshot at a time of maximum shear-stress y at the downstream sensor and (d) base flow.

Figure 7 .

 7 Figure7. Characterization of the attractor without actuation: (a) spectrum of amplitude of the Koopman modes of y, estimated using DFT over a finite timeseries of 300 time units (no windowing), with frequency resolution ∆ω ≈ 2 × 10 -2 ; (b) phase portrait using three x2-velocity probes in the shear-layer at x2 = 0 and x1 = 1/4, 1/2, 3/4 (see figure1); (c) Poincaré section of (b) for u2(1/2, 0) = 0 (crossing from below); the section plane is also indicated in (b). The three plots indicate quasiperiodic behaviour with two incommensurate frequencies (2-torus) ω l ≈ 2.89 and ω h ≈ 11.74.

Figure 8 .Figure 9 .

 89 Figure 8. Spectrogram of the timeseries y(t), using short-time Fourier analysis with Hamming windows of length T = 20 and 50% overlap. In this section, we focus on the unforced flow for t < t1, but further comments regarding the controlled flow are provided in §5.2. Resonant mean-flow modes α, β, γ, δ and b are indicated, according to modal analysis reported in §3.5 and §5.3. The start-up time tm of each controller correction K ′ m is shown with a dashed line.

Figure 10 .

 10 Figure 10. Frequency responses (a) G 0 0 and (b) G b of the mean and base flows, and corresponding sets of poles (c, d) calculated by modal analysis. Eigenvalues are coloured according to criterion (3.7) quantifying the relative contribution of each pole to the input-output relation.The (almost) one-to-one relation between resonance peaks and specific poles is highlighted with dashed vertical lines. These resonant poles can be grouped into a single family for the base flow, highlighted with a grey line in panel (d). This family contains the four unstable modes of the base flow, labelled with greek characters α, β, γ and δ in decreasing order of growth rate and encircled in black. This family of resonant modes can also be identified, by continuity, in the mean flow spectrum. A distinct family of modes a-f , also highlighted with a grey line in panel (c), is responsible for the low-frequency peaks present in |G 0 0 |.

Figure 11 .

 11 Figure 11. Comparison between resolvent-based model G 0 0 (solid red line/solid symbols) and the ROM G 0 0 (dashed black line/open symbols) obtained by subspace identification: frequency response (a) gain and (b) phase (the shaded region corresponds to | arg(G) -arg(G 0 0 )| π); (c) poles α, β, γ, δ.

  |T (iω)| ρ/|W T (iω)|, (4.11) | G 0 0 S(iω)| ρ/|W GS (iω)|, (4.12) |K ′ 1 S(iω)| ρ/|W KS (iω)|,(4.13)

  )

Figure 12 .

 12 Figure 12. Frequency response of the four closed-loop transfer functions of the ROM (black lines): (a) G 0 0 S, (b) S, (c) T , (d) K ′ 1 S. The frequency templates ρ/|WGS|, ρ/|WS|, ρ/|WT | are plotted in blue (ρ/|W K ′ S | is not shown in panel (d) as it is just a constant curve located well above the maximum of |K ′1 S|), with solid lines for ρ = 1 and dashed lines for the optimized-value of the multiplicative constant ρ = 0.6. In panel (b), the shaded region is the forbidden area where the modulus margin constraint (4.20) is violated.

Figure 13 .

 13 Figure13. Timeseries of (a) the perturbation kinetic energy E, the sensor perturbation y and the control signal u = m u ′ m for the controlled flow. The time axis is in log scale to highlight the different time-scales. For each iteration m, the time origin is reset to zero by subtracting tm, the start-up time of the new controller K ′ m , so that the first point is at dt = 4 × 10 -4 in each panel.

Figure 14 .Figure 15 .

 1415 Figure 14. Characterization of dynamical equilibria m = 0, . . . , 5 in the phase plane [∆u2(1/4, 0), ∆u2(3/4, 0)].As the control law is updated, the end-state moves closer and closer to the fixed-point (state 5), which is located at (0, 0). State 0 is a two-torus, states 1-4 are limit cycles (or almost) with decreasing oscillation amplitude and state 5 is the fixed-point. State 2 is in fact an almost degenerate two-torus: the extra dimension of the attractor is so thin that it is hidden by the line width.

5. 3 .

 3 Evolution of the transfer function model G from n to y: gain and poles In figure 16, we follow the evolution of the gain of the transfer function model G characterizing the input-output relation from n to y, during the iteration. All mean flows were computed by averaging over more than 200 time units, in the fully-developped regime. The plant G m m (solid lines), computed using equation (2.32), characterizes the fully-developped flow in feedback loop with the controller K m . At t = t m+1 , a new

Figure 16 .

 16 Figure16. Evolution of the gain |G| during the iterative procedure. Panel (a) shows the initial value of |G| at the beginning of the iteration, at t - 0 . In each subsequent panel, we show the relaxation of G from G m m-1 (dashed lines) at t = t + m to G m m (solid lines) at t = t - m+1 (with t6 → ∞), under the effects of nonlinearity. Same labels as figure15. Vertical lines indicate Koopman frequencies of the permanent state obtained with the controller Km (no unsteady peak for m = 5, because the system is steady). We also indicate which poles of G m m are responsible for each receptivity peak. Mode labels are coloured when their receptivity peak matches with a Koopman frequency of the flow.

Figure 17 .

 17 Figure 17. Evolution of the poles α, β, γ, δ of the ROM transfer function G during the iterative procedure. Solid (respectively open) symbols correspond to poles of G m m (respectively G m+1 m ). Dashed arrows show the downward motion of the poles into the stable half-plane during the shaping of a new closed-loop. Conversely, solid arrows show the upward motion of these poles following nonlinear relaxation of the mean state. Same labels as in figures 15 and 16.

Figure 18 .

 18 Figure 18. Comparison between the gain of the transfer functions G 0 5 = G b associated with the base flow only and G 5 5 , the transfer function associated with the closed-loop between the base flow and the final controller K5.

† Email address for correspondence: colin.leclercq@onera.fr

Acknowledgments

Enlightening discussions with Éric Garnier, Nicolò Fabbiane, Navrose, Olivier Marquet and Yohann Duguet are gratefully acknowledged. The authors are also grateful to Matthew Juniper for sharing his Matlab toolbox for the analysis of nonlinear dynamics (Juniper & Sujith 2018).