2.1 The Navier–Stokes equations as a Lur’e system

We start by considering the non-dimensional, incompressible, Navier–Stokes equations governing the velocity and pressure fields $\boldsymbol{u}$ and $p$

(2.1) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+Re^{-1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$

without external forcing. The Reynolds number $Re:=UL/\unicode[STIX]{x1D708}$ is based on characteristic velocity and length scales $U$ and $L$ of the problem, and on the kinematic viscosity $\unicode[STIX]{x1D708}$ of the fluid. Equations (2.1) are, of course, supplemented by appropriate boundary conditions and an initial condition. Next, we decompose the unknown field $\boldsymbol{q}:=(\boldsymbol{u},p)^{\text{T}}$ ( $^{\text{T}}$ denoting the transpose) into an arbitrary steady incompressible reference field $\boldsymbol{Q}:=(\boldsymbol{U},P)^{\text{T}}$ satisfying the boundary conditions and a perturbation $\boldsymbol{q}^{\prime }:=(\boldsymbol{u}^{\prime },p^{\prime })^{\text{T}}$ with zero boundary conditions, such that $\boldsymbol{q}=\boldsymbol{Q}+\boldsymbol{q}^{\prime }$ . The unsteady perturbation $\boldsymbol{q}^{\prime }$ is governed by the following time-evolution equation

(2.2) $$\begin{eqnarray}{\mathcal{B}}\unicode[STIX]{x2202}_{t}\boldsymbol{q}^{\prime }=[{\mathcal{A}}_{\boldsymbol{ Q}}+\unicode[STIX]{x1D6F9}_{\boldsymbol{Q}}]\boldsymbol{q}^{\prime },\end{eqnarray}$$

where

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{B}}\boldsymbol{q}^{\prime }:=\left(\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{t}\boldsymbol{u}^{\prime }\\ 0\end{array}\right), & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{A}}_{\boldsymbol{Q}}\boldsymbol{q}^{\prime }:=\left(\begin{array}{@{}c@{}}-(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}^{\prime }-(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{U}+Re^{-1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{\prime }-\unicode[STIX]{x1D735}p^{\prime }\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }\end{array}\right), & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F9}_{\boldsymbol{Q}}\boldsymbol{q}^{\prime }:=\left(\begin{array}{@{}c@{}}-(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{U}-\unicode[STIX]{x1D735}P+Re^{-1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{U}-(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}^{\prime }\\ 0\end{array}\right). & \displaystyle\end{eqnarray}$$

The Jacobian ${\mathcal{A}}_{\boldsymbol{Q}}$ and the operator $\unicode[STIX]{x1D6F9}_{\boldsymbol{Q}}$ respectively account for the linear and nonlinear parts of the evolution operator associated with the flow. Since $\boldsymbol{Q}$ is assumed to be time independent, the dynamical system is autonomous; ${\mathcal{A}}_{\boldsymbol{Q}}$ is a linear time-invariant (LTI) operator and $\unicode[STIX]{x1D6F9}_{\boldsymbol{Q}}$ is a static nonlinearity. Taking the Laplace transform of (2.2) leads to the relation

(2.6) $$\begin{eqnarray}\boldsymbol{q}^{\prime }(s)={\mathcal{R}}_{\boldsymbol{ Q}}(s)[\unicode[STIX]{x1D6F9}_{\boldsymbol{Q}}\boldsymbol{q}^{\prime }(s)+{\mathcal{B}}\boldsymbol{q}_{t=0}^{\prime }],\end{eqnarray}$$

where $s$ is the Laplace variable, $\boldsymbol{q}_{t=0}^{\prime }$ is the initial condition on the perturbation and

(2.7) $$\begin{eqnarray}{\mathcal{R}}_{\boldsymbol{Q}}(s):=(s{\mathcal{B}}-{\mathcal{A}}_{\boldsymbol{Q}})^{-1}\end{eqnarray}$$

is called the resolvent operator about the reference field $\boldsymbol{Q}$ . Such interconnection of an LTI operator ${\mathcal{R}}_{\boldsymbol{Q}}$ with a static nonlinearity $\unicode[STIX]{x1D6F9}_{\boldsymbol{Q}}$ , illustrated in figure 2, is called a Lur’e system in nonlinear control theory (Khalil Reference Khalil2002; Sharma Reference Sharma2009; McKeon & Sharma Reference McKeon and Sharma2010).

Figure 2.

The Navier–Stokes equations as a Lur’e system.

The resolvent operator is defined for any $s$ which is not an eigenvalue of the Jacobian operator, i.e. a solution of

(2.8) $$\begin{eqnarray}(s{\mathcal{B}}-{\mathcal{A}}_{\boldsymbol{Q}})\boldsymbol{q}^{\prime }=0,\end{eqnarray}$$

with non-zero $\boldsymbol{q}^{\prime }$ . Interestingly, the output $\boldsymbol{q}^{\prime }$ remains bounded even if ${\mathcal{A}}_{\boldsymbol{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 ${\mathcal{A}}_{\boldsymbol{Q}}$ has no marginally unstable eigenvalues. The choice of the arbitrary reference field $\boldsymbol{Q}$ is discussed in the next paragraph.

2.2 Open-loop and closed-loop transfer functions of fully developed flows

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 $G_{\boldsymbol{Q}}$ from $n$ to $y$ and an output disturbance $w$ in (b) the open loop and (d) the closed loop.

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

(2.9) $$\begin{eqnarray}y(t):=\langle \boldsymbol{C},\boldsymbol{q}^{\prime }(t)\rangle +Y,\end{eqnarray}$$

where

(2.10) $$\begin{eqnarray}\langle \boldsymbol{q}_{1},\boldsymbol{q}_{2}\rangle :=\int _{{\mathcal{D}}}\boldsymbol{q}_{1}^{\ast }\boldsymbol{\cdot }\boldsymbol{q}_{2}\,\text{d}{\mathcal{D}}\end{eqnarray}$$

defines an inner product over the domain ${\mathcal{D}}$ and $^{\ast }$ denotes complex conjugation. The output $y$ is linear with respect to the fluctuation $\boldsymbol{q}^{\prime }$ 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

(2.11) $$\begin{eqnarray}u:=Ky\end{eqnarray}$$

is due to an LTI controller of transfer function $K(s)$ . With external forcing, the output $\boldsymbol{q}^{\prime }$ is respectively given by

(2.12) $$\begin{eqnarray}\boldsymbol{q}^{\prime }={\mathcal{R}}_{\boldsymbol{ Q}}[\unicode[STIX]{x1D6F9}_{\boldsymbol{Q}}\boldsymbol{q}^{\prime }+{\mathcal{B}}\boldsymbol{q}_{t=0}^{\prime }+\boldsymbol{B}n]\end{eqnarray}$$

in the open-loop case, and

(2.13) $$\begin{eqnarray}\boldsymbol{q}^{\prime }=({\mathcal{I}}-{\mathcal{R}}_{\boldsymbol{ Q}}\boldsymbol{B}K\langle \boldsymbol{C},\boldsymbol{\cdot }\rangle )^{-1}{\mathcal{R}}_{\boldsymbol{Q}}\left[\unicode[STIX]{x1D6F9}_{\boldsymbol{Q}}\boldsymbol{q}^{\prime }+{\mathcal{B}}\boldsymbol{q}_{t=0}^{\prime }+\boldsymbol{B}n+\boldsymbol{B}K{\displaystyle \frac{Y}{s}}\right]\end{eqnarray}$$

in the closed-loop case. The latter expression can be rewritten as

(2.14) $$\begin{eqnarray}\boldsymbol{q}^{\prime }={\mathcal{R}}_{\boldsymbol{ Q}}^{K}\left[\unicode[STIX]{x1D6F9}_{\boldsymbol{Q}}\boldsymbol{q}^{\prime }+{\mathcal{B}}\boldsymbol{q}_{t=0}^{\prime }+\boldsymbol{B}n+\boldsymbol{B}K{\displaystyle \frac{Y}{s}}\right],\end{eqnarray}$$

where

(2.15) $$\begin{eqnarray}{\mathcal{R}}_{\boldsymbol{Q}}^{K}(s):=(s{\mathcal{B}}-{\mathcal{A}}_{\boldsymbol{ Q}}^{K}(s))^{-1}\end{eqnarray}$$

is the resolvent operator associated with the Jacobian operator

(2.16) $$\begin{eqnarray}{\mathcal{A}}_{\boldsymbol{Q}}^{K}(s):={\mathcal{A}}_{\boldsymbol{ Q}}+K(s)\boldsymbol{B}\langle \boldsymbol{C},.\rangle\end{eqnarray}$$

of the coupled system fluid/controller. Because of the nonlinear feedback term, the concept of a transfer function from $n$ to $\boldsymbol{q}^{\prime }$ is ill defined. However, using the resolvent operators introduced above and lumping all the terms independent of $n$ in an output disturbance term $\boldsymbol{w}$ , we can recast the governing equations (2.12) and (2.14) in a form that resembles input–output relations:

(2.17a,b ) $$\begin{eqnarray}\boldsymbol{q}^{\prime }={\mathcal{R}}_{\boldsymbol{ Q}}\boldsymbol{B}n+\boldsymbol{w}\quad \text{and}\quad \boldsymbol{q}^{\prime }={\mathcal{R}}_{\boldsymbol{ Q}}^{K}\boldsymbol{B}n+\boldsymbol{w}\end{eqnarray}$$

in the open-loop and closed-loop cases respectively. Similar to the unforced case, the output $\boldsymbol{q}^{\prime }$ is bounded despite the potential presence of unstable eigenvalues in ${\mathcal{A}}_{\boldsymbol{Q}}$ or ${\mathcal{A}}_{\boldsymbol{Q}}^{K}$ , thanks to the nonlinear feedback term (included in $\boldsymbol{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 $\boldsymbol{Q}$ . The passivity theorem states that if the two interconnected blocks ${\mathcal{R}}_{\boldsymbol{Q}}^{K}$ and $\unicode[STIX]{x1D6F9}_{\boldsymbol{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 $\boldsymbol{Q}=\boldsymbol{q}_{b}$ , the static nonlinearity

(2.18) $$\begin{eqnarray}\unicode[STIX]{x1D6F9}_{\boldsymbol{q}_{b}}\boldsymbol{q}^{\prime }=\left(\begin{array}{@{}c@{}}-(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}^{\prime }\\ 0\end{array}\right)\end{eqnarray}$$

happens to be conservative, hence passive. Therefore, a possible strategy to reach the base flow from the regime of fully developed 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 ${\mathcal{A}}_{\boldsymbol{q}_{b}}$ . This strategy has been pursued in recent works by Sharma et al. (Reference Sharma, Morrison, McKeon, Limebeer, Koberg and Sherwin2011) and Heins et al. (Reference Heins, Jones and Sharma2016), 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 developed regime. McKeon & Sharma (Reference McKeon and Sharma2010) therefore proposed to consider the mean flow $\overline{\boldsymbol{q}}$ as an appropriate alternative choice of reference field $\boldsymbol{Q}$ in (2.6), in order to define the unforced Navier–Stokes as a Lur’e system. In that case, the nonlinear output

(2.19) $$\begin{eqnarray}\unicode[STIX]{x1D6F9}_{\overline{\boldsymbol{q}}}\boldsymbol{q}^{\prime }=\left(\begin{array}{@{}c@{}}\overline{(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}^{\prime }}-(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}^{\prime }\\ 0\end{array}\right)\end{eqnarray}$$

was interpreted as an endogenous forcing term. We decide to follow the same approach here and also choose

(2.20) $$\begin{eqnarray}\boldsymbol{Q}:=\overline{\boldsymbol{q}}\end{eqnarray}$$

to model the flow receptivity to the exogenous forcing field $\boldsymbol{B}$ .

In the case of unforced flows, the mean flow is unambiguously defined as the temporal average of $\boldsymbol{q}$ in the fully developed 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: (i) that associated with the initial attractor characterizing the unforced system at $t=0$ , and (ii) that characterizing the final attractor of the forced system as $t\rightarrow \infty$ . During the transient, the mean flow is not defined. If we want to extend the framework of McKeon & Sharma (Reference McKeon and Sharma2010) 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

(2.21) $$\begin{eqnarray}y=G_{\overline{\boldsymbol{q}}}n+w,\end{eqnarray}$$

where

(2.22) $$\begin{eqnarray}G_{\overline{\boldsymbol{q}}}(s):=\langle \boldsymbol{C},{\mathcal{R}}_{\overline{\boldsymbol{q}}}(s)\boldsymbol{B}\rangle\end{eqnarray}$$

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

(2.23) $$\begin{eqnarray}G_{\overline{\boldsymbol{q}}}^{K}(s):=\langle \boldsymbol{C},{\mathcal{R}}_{\overline{\boldsymbol{q}}}^{K}(s)\boldsymbol{B}\rangle\end{eqnarray}$$

and depends on both $\overline{\boldsymbol{q}}$ and the controller $K$ . This pseudo-transfer function can also be computed from the expression

(2.24) $$\begin{eqnarray}G_{\overline{\boldsymbol{q}}}^{K}={\displaystyle \frac{G_{\overline{\boldsymbol{q}}}}{1-G_{\overline{\boldsymbol{q}}}K}},\end{eqnarray}$$

which corresponds to the linear plant $G_{\overline{\boldsymbol{q}}}$ in the feedback loop with the controller $K$ , as can be seen in figure 3(d). The reference field $\overline{\boldsymbol{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.

2.3 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}^{\prime }$ such that the total controller $K_{2}:=K_{1}+K_{2}^{\prime }$ is able to damp the oscillations further? In general, from a sequence of $m$ successive controller corrections

(2.25) $$\begin{eqnarray}K_{m}:=\left\{\begin{array}{@{}ll@{}}0,\quad & \text{if }m=0,\\ K_{1}^{\prime }+K_{2}^{\prime }+\cdots +K_{m}^{\prime },\quad & \text{if }m>0,\end{array}\right.\end{eqnarray}$$

leading to a permanent oscillatory regime, how can we design the next controller correction $K_{m+1}^{\prime }$ such that $K_{m+1}=K_{m}+K_{m+1}^{\prime }$ leads to even smaller oscillation amplitude?

The key idea is to consider an input–output model based on the mean flow $\overline{\boldsymbol{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 ${\mathcal{R}}_{\overline{\boldsymbol{q}}_{m}}^{K_{m}}$ , involving both $\overline{\boldsymbol{q}}_{m}$ and $K_{m}$ . The corresponding transfer function from $n$ to $y$ is then given by $G_{\overline{\boldsymbol{q}}_{m}}^{K_{m}}$ . To simplify notations, we denote from now on

(2.26) $$\begin{eqnarray}G_{i}^{j}:=G_{\overline{\boldsymbol{q}}_{i}}^{K_{j}}.\end{eqnarray}$$

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}^{\prime }$ such that the target plant

(2.27) $$\begin{eqnarray}G_{m}^{m+1}={\displaystyle \frac{G_{m}^{m}}{1-G_{m}^{m}K_{m+1}^{\prime }}}\end{eqnarray}$$

has some desired properties, leading to smaller oscillations of the output $y$ . When the additional controller $K_{m+1}^{\prime }$ is added in the feedback loop, the flow transiently evolves until a new dynamical equilibrium is reached, characterized by a new mean state $\overline{\boldsymbol{q}}_{m+1}$ . Since the effective controller is now $K_{m+1}=K_{m}+K_{m+1}^{\prime }$ , 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.

Figure 4.

Evolution of the block diagram representation of the system during the iterative procedure: (a) $t=t_{1}^{+}$ , $K_{1}^{\prime }$ has just 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 $K_{1}:=K_{1}^{\prime }$ has converged to a new mean state; (c) $t=t_{2}^{+}$ , $K_{2}^{\prime }$ has just been added so the flow alone is still characterized by the plant $G_{1}^{0}$ ; (d) $t=t_{2}^{-}$ the feedback loop between the fluid and the controller $K_{2}:=K_{1}^{\prime }+K_{2}^{\prime }$ has converged to a new mean state.

The momentum equation with feedback reads

(2.28) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+Re^{-1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\boldsymbol{B}_{\boldsymbol{ u}}\int _{0}^{t}\unicode[STIX]{x1D705}_{m}(\unicode[STIX]{x1D70F})y(t-\unicode[STIX]{x1D70F})\,\text{d}\unicode[STIX]{x1D70F},\end{eqnarray}$$

where $\boldsymbol{B}_{\boldsymbol{u}}$ is the component of $\boldsymbol{B}$ forcing the momentum equation and $\unicode[STIX]{x1D705}_{m}$ is the impulse response of the controller $K_{m}$ . If a fixed point $\boldsymbol{q}_{\infty }:=(\boldsymbol{u}_{\infty },p_{\infty })$ is reached in the limit $t\rightarrow \infty$ , then it must satisfy

(2.29) $$\begin{eqnarray}(\boldsymbol{u}_{\infty }\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}_{\infty }=-\unicode[STIX]{x1D735}p_{\infty }+Re^{-1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{\infty }+\boldsymbol{B}_{\boldsymbol{u}}K_{m}(0)(\langle \boldsymbol{C},\boldsymbol{q}_{\infty }-\boldsymbol{Q}\rangle +Y).\end{eqnarray}$$

In order to reach the base flow, i.e.  $\boldsymbol{q}_{\infty }=\boldsymbol{q}_{b}$ , we must ensure that $\boldsymbol{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

(2.30) $$\begin{eqnarray}K_{j}^{\prime }(0)=0\end{eqnarray}$$

at every step $j$ . In the special case where the sensor measures a deviation with respect to the base flow by satisfying

(2.31) $$\begin{eqnarray}Y=\langle \boldsymbol{C},\boldsymbol{Q}-\boldsymbol{q}_{b}\rangle ,\end{eqnarray}$$

then $\boldsymbol{q}_{b}$ remains a fixed point of (2.29) even if $K_{m}(0)\neq 0$ , but there may be alternative fixed points $\boldsymbol{q}_{\infty }$ as well. It is therefore crucial to ‘disconnect’ the controller in the steady regime by imposing (2.30); then $\boldsymbol{q}_{b}$ is for sure the only fixed point for any choice of sensor $(\boldsymbol{C},Y)$ .

Although the iterative method bears resemblance to the gain scheduling approach of Khalil (Reference Khalil2002) and Högberg, Bewley & Henningson (Reference Högberg, Bewley and Henningson2003), there are fundamental differences which are emphasized later in the discussion, see § 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 (i) the system is in a dynamical equilibrium characterized by $\overline{\boldsymbol{q}}_{m}$ and $K_{m}$ , or (ii) the system is out of equilibrium because a controller correction $K_{m+1}^{\prime }$ has just been added, but the flow has not yet evolved so the mean flow $\overline{\boldsymbol{q}}_{m}$ remains relevant over some period of time. Denoting $t_{m+1}$ the instant when the controller correction $K_{m+1}^{\prime }$ is added, we have $G=G_{m}^{m}$ at $t=t_{m+1}^{-}$ and $G=G_{m}^{m+1}$ 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 (2.24) and (2.25), each transfer function $G_{i}^{j}$ can be seen as the interconnection

(2.32) $$\begin{eqnarray}G_{i}^{j}={\displaystyle \frac{G_{i}^{0}}{1-G_{i}^{0}(K_{1}^{\prime }+\cdots +K_{j}^{\prime })}}\end{eqnarray}$$

between the plant $G_{i}^{0}$ characterizing the mean flow alone (even though it cannot be sustained without controllers for $i\neq 0$ ), and the controller corrections $K_{1}^{\prime }$ to $K_{j}^{\prime }$ plugged in parallel. The block diagrams at times $t=t_{1}^{+}$ , $t_{2}^{\pm }$ and $t_{3}^{-}$ are shown in figure 4.

2.4 Relation between poles of the transfer operators $G_{i}^{j}$ and instability modes

In the present context, modal analysis of the Jacobian operator ${\mathcal{A}}_{i}^{j}$ is used to find the pole singularities of the resolvent operator ${\mathcal{R}}_{i}^{j}$ (with obvious notations), and therefore of the transfer function $G_{i}^{j}$ . We are concerned with the solutions of the generalized eigenvalue problem

(2.33) $$\begin{eqnarray}(s_{l}{\mathcal{B}}-{\mathcal{A}}_{i}^{j}(s_{l}))\boldsymbol{q}_{l}^{\prime }=0,\end{eqnarray}$$

with complex eigenvalues $s_{l}=\text{i}\unicode[STIX]{x1D714}_{l}+\unicode[STIX]{x1D70E}_{l}$ of growth rate $\unicode[STIX]{x1D70E}_{l}$ , frequency $\unicode[STIX]{x1D714}_{l}$ and eigenmodes $\boldsymbol{q}_{l}^{\prime }$ (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 $\overline{\boldsymbol{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 ${\mathcal{A}}_{i}^{j}$ may be seen as a perturbation of operator ${\mathcal{A}}_{i}^{0}$ with a correction term proportional to $K_{j}$ according to (2.16), i.e.

(2.34) $$\begin{eqnarray}{\mathcal{A}}_{i}^{j}={\mathcal{A}}_{i}^{0}+K_{j}\boldsymbol{B}\langle \boldsymbol{C},\boldsymbol{\cdot }\rangle .\end{eqnarray}$$

A continuous connection can therefore be established between the eigenvalues of the two operators. The Jacobian operator ${\mathcal{A}}_{i}^{0}$ may in turn be seen as a perturbation of the Jacobian operator about the base flow ${\mathcal{A}}_{\boldsymbol{q}_{b}}^{0}$ , with a correction term due to mean-flow distortion $\unicode[STIX]{x0394}\overline{\boldsymbol{q}}_{i}$ , where

(2.35) $$\begin{eqnarray}\unicode[STIX]{x0394}\boldsymbol{q}:=\boldsymbol{q}-\boldsymbol{q}_{b}.\end{eqnarray}$$

Therefore, all instability modes of the base flow may be unambiguously identified with some poles of $G_{i}^{j}$ , thereby providing the latter with a physical meaning. This key advantage of our resolvent-based model will be illustrated in § 3.5.

3.1 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 (Barbagallo et al. Reference Barbagallo, Sipp and Schmid2009, Reference Barbagallo, Sipp and Schmid2011; Dergham et al. Reference Dergham, Sipp, Robinet and Barbagallo2011; Poussot-Vassal & Sipp Reference Poussot-Vassal and Sipp2015; Schmid & Sipp Reference Schmid and Sipp2016; Sipp & Schmid Reference Sipp and Schmid2016) following Sipp & Lebedev (Reference Sipp and Lebedev2007). The reader is referred to this initial paper for a comprehensive description of the numerical set-up. 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 $\boldsymbol{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 $\unicode[STIX]{x2202}_{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 ${\sim}200\,000$ $(P1b,P1b,P1)$ Taylor–Hood elements, yielding a total of ${\sim}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 (Reference Sipp and Schmid2016), 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 $\boldsymbol{q}_{b}$ – which was computed using Newton’s iteration (Sipp & Schmid Reference Sipp and Schmid2016) – 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 (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1997), coupled with the multifrontal sparse LU solver MUMPS (Amestoy et al. Reference Amestoy, Duff, Koster and L’Excellent2001, Reference Amestoy, Guermouche, L’Excellent and Pralet2006). The shift-invert mode was used to find eigenvalues of the Jacobian operator in the vicinity of a given shift.

3.2 Actuator and flow measurements

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

(3.1) $$\begin{eqnarray}\boldsymbol{B}=\left[0,\unicode[STIX]{x1D702}\exp \left(-{\displaystyle \frac{(x_{1}-x_{1}^{0})^{2}+(x_{2}-x_{2}^{0})^{2}}{2\unicode[STIX]{x1D70E}^{2}}}\right),0\right]^{\text{T}},\end{eqnarray}$$

centred around $(x_{1}^{0},x_{2}^{0})=(-0.1,0.02)$ , with a scalar $\unicode[STIX]{x1D702}$ such that $\langle \boldsymbol{B},\boldsymbol{B}\rangle =1$ . With the choice of spatial variance $\unicode[STIX]{x1D70E}=0.0849$ , the volume force reaches $50\,\%$ of its maximum value at a radial distance of ${\approx}0.1$ from $(x_{1}^{0},x_{2}^{0})$ , and $1\,\%$ at ${\approx}0.26$ .

In order to characterize the strength of the nonlinearity, we define different measures of the distance between the instantaneous field $\boldsymbol{q}$ and the base flow $\boldsymbol{q}_{b}$ , based on the perturbation field $\unicode[STIX]{x0394}\boldsymbol{q}$ . The downstream-edge sensor is modelled by a pair $(\boldsymbol{C},Y)$ such that

(3.2) $$\begin{eqnarray}y=\left.\int _{x_{1}=1}^{1.1}\unicode[STIX]{x2202}_{2}\unicode[STIX]{x0394}u_{1}\right|_{x_{2}=0}\,\text{d}x_{1}\end{eqnarray}$$

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}^{\prime }(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

(3.3) $$\begin{eqnarray}E:={\displaystyle \frac{\langle \unicode[STIX]{x0394}\boldsymbol{q},{\mathcal{B}}\unicode[STIX]{x0394}\boldsymbol{q}\rangle }{\langle \boldsymbol{q}_{b},{\mathcal{B}}\boldsymbol{q}_{b}\rangle }}.\end{eqnarray}$$

3.3 Transient dynamics

Figure 5.

Time series of (a–c) perturbation kinetic energy $E$ and (d–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 5 shows time series of the global and local measures $E$ and $y$ of the uncontrolled flow. We recall that the simulation is initialized with the base flow plus 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 $\overline{E}$ is of the order of $2\,\%$ (see figure 5 c), which corresponds to perturbation velocities of the order of $0.1$ .

3.4 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.

Comparison between (a,c) state 0 and (b,d) the base flow. (a,b) Velocity 2-norm of (a) uncontrolled mean flow $\overline{\boldsymbol{u}}_{0}$ (computed by averaging over 200 time units, after the transient) and (b) base flow $\boldsymbol{u}_{b}$ , with associated contours of the streamfunction $\unicode[STIX]{x1D713}$ . The streamfunction is defined by the relations $\unicode[STIX]{x2202}_{2}\unicode[STIX]{x1D713}=v_{1}$ , $\unicode[STIX]{x2202}_{1}\unicode[STIX]{x1D713}=-v_{2}$ and $\unicode[STIX]{x1D713}=0$ on the lower wall (solid contours for positive increments $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}=0.05$ from $\unicode[STIX]{x1D713}=0$ ; dashed contours for negative increments $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}=-0.01$ from $\unicode[STIX]{x1D713}=0$ ). (c,d) Vorticity fields: (c) snapshot at a time of maximum shear stress $y$ at the downstream sensor and (d) base flow.

Figure 6(c) shows a vorticity snapshot of state 0 yielding maximal instantaneous shear stress $y$ , compared to the steady vorticity field of the base flow in (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 $\overline{\langle \boldsymbol{q}^{\prime },\boldsymbol{q}^{\prime }\rangle }/\langle \overline{\boldsymbol{q}}_{0},\overline{\boldsymbol{q}}_{0}\rangle \approx 1.1\times 10^{-3}$ is very low, so the mean flow can be considered as a good approximation of the full flow at any instant.

Figure 7.

Characterization of the attractor without actuation: (a) spectrum of amplitude of the Koopman modes of $y$ , estimated using DFT over a finite time series of $300$ time units (no windowing), with frequency resolution $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}\approx 2\times 10^{-2}$ ; (b) phase portrait using three $x_{2}$ -velocity probes in the shear layer at $x_{2}=0$ and $x_{1}=1/4,1/2,3/4$ (see figure 1); (c) Poincaré section of (b) for $u_{2}(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) $\unicode[STIX]{x1D714}_{l}\approx 2.89$ and $\unicode[STIX]{x1D714}_{h}\approx 11.74$ .

The dynamical state may be further characterized in the frequency domain by considering the spectrum of Koopman modes. For fully developed flows, these can be computed using harmonic averages

(3.4) $$\begin{eqnarray}\hat{\boldsymbol{q}}^{\prime }(\unicode[STIX]{x1D714}):=\lim _{T\rightarrow \infty }{\displaystyle \frac{1}{T}}\int _{0}^{T}\boldsymbol{q}^{\prime }(t)\text{e}^{-\text{i}\unicode[STIX]{x1D714}t}\,\text{d}t\end{eqnarray}$$

(Arbabi & Mezić Reference Arbabi and Mezić2017). The limit is non-zero for Koopman eigenfrequencies only, and yields the associated Koopman modes of the observable $\boldsymbol{q}^{\prime }$ . In practice, the spectrum of amplitude of the Koopman modes may be approximated using discrete Fourier transform (DFT) on a finite time series 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 three-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 $\unicode[STIX]{x0394}u_{2}(1/2,0)=0$ from below yields a Poincaré map in the $[\unicode[STIX]{x0394}u_{2}(1/4,0),\unicode[STIX]{x0394}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 $\unicode[STIX]{x1D714}_{l}$ and $\unicode[STIX]{x1D714}_{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:

(3.5) $$\begin{eqnarray}\{\unicode[STIX]{x1D714}_{\boldsymbol{k}}\}:=\{k_{l}\unicode[STIX]{x1D714}_{l}+k_{h}\unicode[STIX]{x1D714}_{h}\mid \boldsymbol{k}:=(k_{l},k_{h})\in \mathbb{Z}^{2}\}.\end{eqnarray}$$

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.

Figure 8.

Spectrogram of the time series $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<t_{1}$ , but further comments regarding the controlled flow are provided in § 5.2. Resonant mean-flow modes $\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FF}$ and $b$ are indicated, according to modal analysis reported in §§ 3.5 and 5.3. The start-up time $t_{m}$ of each controller correction $K_{m}^{\prime }$ is shown with a dashed line.

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 $\unicode[STIX]{x1D714}=10.9$ (Barbagallo et al. Reference Barbagallo, Sipp and Schmid2009), 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 $\unicode[STIX]{x1D714}_{h}$ associated with shear-layer oscillations. During the transient, we also observe the development of a low-frequency peak at $\unicode[STIX]{x1D714}_{l}\approx 2.89$ . The flow structures associated with $\unicode[STIX]{x1D714}_{l}$ and $\unicode[STIX]{x1D714}_{h}$ are shown in figure 9(a,b). Panel (b) confirms that the Koopman mode at $\unicode[STIX]{x1D714}_{h}$ is associated with a shear-layer mode. On the other hand, the low-frequency mode in (a) is localized within the cavity. The frequency is low because the recirculating mean flow within the cavity is slow. The two ‘fundamental’ frequencies $\unicode[STIX]{x1D714}_{h}$ and $\unicode[STIX]{x1D714}_{l}$ give rise to harmonics and interaction peaks visible in the spectrogram and the spectrum of state 0 in figure 7(a).

Figure 9.

Comparison between real part of streamwise velocity component of Koopman modes at (a) $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{l}$ , (b) $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{h}$ and eigenmodes (c) $b$ , (d) $\unicode[STIX]{x1D6FC}$ of the uncontrolled mean flow. Koopman modes $\hat{\boldsymbol{q}}^{\prime }$ were approximated using DFT on a time series of 4013 snapshots equally sampled over a total duration of $128.416$ time units. The complex amplitude of the mean-flow eigenmodes was tuned to minimize the $L_{2}$ -error with respect to the Koopman modes. The colour map spans from the minimum to the maximum value in each plot, so bounds can be compared between Koopman modes and eigenmodes.

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 (d). This family contains the four unstable modes of the base flow, labelled with Greek characters $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ , $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FF}$ 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 (c), is responsible for the low-frequency peaks present in $|G_{0}^{0}|$ .

3.5 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 $\overline{\boldsymbol{q}}_{0}$ correspond to poles of the transfer function $G_{0}^{0}$ , which can be related to instability modes of the base flow $\boldsymbol{q}_{b}$ . On the top panels, we plot the gain of the frequency response associated with the unforced mean flow $G_{0}^{0}(\text{i}\unicode[STIX]{x1D714})$ (a) and base flow $G_{b}(\text{i}\unicode[STIX]{x1D714})$ (b), the latter being defined as

(3.6) $$\begin{eqnarray}G_{b}:=\langle \boldsymbol{C},{\mathcal{R}}_{\boldsymbol{q}_{b}}\boldsymbol{B}\rangle .\end{eqnarray}$$

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 (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 mean-flow or base-flow eigenvalues. We therefore plot the spectra of the associated Jacobian operators in (c,d), below the frequency responses. In order to identify which eigenvalues are responsible for the receptivity peaks, we colour the spectra by the criterion

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{j}={\displaystyle \frac{|\langle \boldsymbol{C},\boldsymbol{q}_{j}^{\prime }\rangle ||\langle \boldsymbol{p}_{j}^{\prime },\boldsymbol{B}\rangle |}{|\unicode[STIX]{x1D70E}_{j}|}},\end{eqnarray}$$

quantifying the contribution of each mode to the input–output relation (Antoulas Reference Antoulas2005; Bagheri et al. Reference Bagheri, Brandt and Henningson2009; Barbagallo et al. Reference Barbagallo, Sipp and Schmid2009). The criterion takes into account the distance $1/|\unicode[STIX]{x1D70E}_{l}|$ of each eigenvalue to the imaginary axis, the observability of the direct mode $|\langle \boldsymbol{C},\boldsymbol{q}_{l}^{\prime }\rangle |$ and the controllability of the adjoint mode $|\langle \boldsymbol{p}_{l}^{\prime },\boldsymbol{B}\rangle |$ . We recall that the adjoint modes $(\boldsymbol{p}_{l}^{\prime })$ are solutions of the eigenvalue problem $(\unicode[STIX]{x1D707}_{l}{\mathcal{B}}-\widetilde{{\mathcal{A}}}_{i}^{j}(\unicode[STIX]{x1D707}_{l}))\boldsymbol{p}_{l}^{\prime }=0$ , where $\widetilde{{\mathcal{A}}}_{i}^{j}$ is the adjoint of ${\mathcal{A}}_{i}^{j}$ 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.  $\unicode[STIX]{x1D707}_{l}=s_{l}^{\ast }$ , and adjoint modes satisfy the biorthogonality relation $\langle \boldsymbol{p}_{k}^{\prime },\boldsymbol{q}_{l}^{\prime }\rangle =\unicode[STIX]{x1D6FF}_{kl}$ with the normalized direct modes $(\boldsymbol{q}_{l}^{\prime })$ . 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 (c,d), intended to guide the eye. Four poles, labelled $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ , $\unicode[STIX]{x1D6FE}$ , $\unicode[STIX]{x1D6FF}$ are responsible for the peaks at $\unicode[STIX]{x1D714}>5$ in $G_{0}^{0}$ , corresponding to oscillations at $\unicode[STIX]{x1D714}_{h}$ , $\unicode[STIX]{x1D714}_{h}\pm \unicode[STIX]{x1D714}_{l}$ and $\unicode[STIX]{x1D714}_{h}+2\unicode[STIX]{x1D714}_{l}$ . The other family of mean-flow modes $a,b,c,d,e,f$ are responsible for the low-frequency peaks at $\unicode[STIX]{x1D714}<5$ , among which the peak at $\unicode[STIX]{x1D714}_{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 $\unicode[STIX]{x1D714}_{h}$ is connected to the most unstable mode of the base flow. Therefore, mode $\unicode[STIX]{x1D6FC}$ is the mean-flow analogue of that specific base-flow eigenmode. The other mean-flow modes $\unicode[STIX]{x1D6FD}$ , $\unicode[STIX]{x1D6FE}$ , $\unicode[STIX]{x1D6FF}$ 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 $\unicode[STIX]{x1D714}_{l}$ in $G_{0}^{0}$ does not have an obvious base-flow analogue. In fact, the emergence of mode $b$ in the mean-flow 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 $\{\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FF}\}$ is a multiple of the frequency $\unicode[STIX]{x1D714}_{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 (Cherubini, Robinet & De Palma Reference Cherubini, Robinet and De Palma2010).

In figure 9, we compare the structure of the mean-flow eigenmodes $\unicode[STIX]{x1D6FC}$ and $b$ (c,d) with the Koopman modes of the unforced flow at $\unicode[STIX]{x1D714}_{h}$ and $\unicode[STIX]{x1D714}_{l}$ (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 (Turton et al. Reference Turton, Tuckerman and Barkley2015; Mantič-Lugo & Gallaire Reference Mantič-Lugo and Gallaire2016), the so-called RZIF property (real-zero imaginary-frequency) in Turton et al. (Reference Turton, Tuckerman and Barkley2015). We highlight the fact that the RZIF property is not verified here: modes $b$ , $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ , $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FF}$ have non-zero, positive growth rates, yet they are associated with the correct frequencies and Koopman modes. This was already noticed by Mettot, Renac & Sipp (Reference Mettot, Renac and Sipp2014a ) and Meliga (Reference Meliga2017) in open-cavity 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 (Mezić Reference Mezić2013), as previously said in § 3.4.

Finally, we observe in figure 10 receptivity peaks at low frequencies different from $\unicode[STIX]{x1D714}_{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

(3.8) $$\begin{eqnarray}\boldsymbol{f}^{\prime }:=\unicode[STIX]{x1D6F9}_{\overline{\boldsymbol{q}}}\boldsymbol{q}^{\prime }\end{eqnarray}$$

due to the nonlinear feedback operator. Since the operator $\unicode[STIX]{x1D6F9}_{\overline{\boldsymbol{q}}}$ is time invariant, its output $\boldsymbol{f}^{\prime }$ is quasiperiodic with the same discrete set of frequencies as its input $\boldsymbol{q}^{\prime }$ , i.e.  $\{\unicode[STIX]{x1D714}_{\boldsymbol{k}}\}$ defined in (3.5). More precisely, we can introduce the sets of Koopman modes

(3.9a,b ) $$\begin{eqnarray}\{\hat{\boldsymbol{q}}_{\boldsymbol{k}}^{\prime }\}:=\{\hat{\boldsymbol{q}}^{\prime }(\unicode[STIX]{x1D714}_{\boldsymbol{ k}})\},\quad \{\hat{\boldsymbol{f}}_{\boldsymbol{k}}^{\prime }\}:=\{\hat{\boldsymbol{f}}^{\prime }(\unicode[STIX]{x1D714}_{\boldsymbol{ k}})\}\end{eqnarray}$$

respectively associated with $\boldsymbol{q}^{\prime }$ and $\boldsymbol{f}^{\prime }$ . In the fully developed regime, the governing equations (2.2) can be recast in the form of a harmonic balance relation (Khalil Reference Khalil2002)

(3.10) $$\begin{eqnarray}\hat{\boldsymbol{q}}_{\boldsymbol{k}}^{\prime }={\mathcal{R}}_{\overline{\boldsymbol{q}}}(\text{i}\unicode[STIX]{x1D714}_{\boldsymbol{ k}})\hat{\boldsymbol{f}}_{\boldsymbol{k}}^{\prime },\end{eqnarray}$$

between these two sets of Koopman modes (harmonic averages are zero otherwise). Clearly, strong receptivity at frequencies outside the set $\{\unicode[STIX]{x1D714}_{\boldsymbol{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.

4.1 Model reduction

We convert the frequency response $G_{0}^{0}$ in the form of a state-space model $(\widetilde{\boldsymbol{A}},\widetilde{\boldsymbol{B}},\widetilde{\boldsymbol{C}},\widetilde{\boldsymbol{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.  $\widetilde{\boldsymbol{D}}=0$ ) reduced-order model (ROM) of order 16 is obtained using subspace identification methods (Liu, Jacques & Miller Reference Liu, Jacques and Miller1996; McKelvey, Akcay & Ljung Reference McKelvey, Akcay and Ljung1996). For subsequent iterations, the order of the ROM varies between 16 and 32. The frequency response of the ROM is given by

(4.1) $$\begin{eqnarray}\widetilde{G}_{0}^{0}(\text{i}\unicode[STIX]{x1D714})=\widetilde{\boldsymbol{C}}(\text{i}\unicode[STIX]{x1D714}\boldsymbol{I}-\widetilde{\boldsymbol{A}})^{-1}\widetilde{\boldsymbol{B}}\end{eqnarray}$$

and the associated poles correspond to eigenvalues of $\widetilde{\boldsymbol{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 $\widetilde{G}_{0}^{0}(\text{i}\unicode[STIX]{x1D714})$ (black dotted curve) coincides with the original model $G_{0}^{0}$ (solid red curve) over the frequency range $5<\unicode[STIX]{x1D714}<20$ . This range contains the main resonance peaks at $\unicode[STIX]{x1D714}_{h}$ , $\unicode[STIX]{x1D714}_{h}\pm \unicode[STIX]{x1D714}_{l}$ and $\unicode[STIX]{x1D714}_{h}+2\unicode[STIX]{x1D714}_{l}$ , but not the low-frequency ones at $\unicode[STIX]{x1D714}<5$ which we choose to ignore in order to ensure sufficient order reduction. In (c), we check that the poles $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ , $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FF}$ corresponding to the main resonance peaks are correctly approximated by some eigenvalues of $\widetilde{\boldsymbol{A}}$ .

Figure 11.

Comparison between resolvent-based model $G_{0}^{0}$ (solid red line/solid symbols) and the ROM $\widetilde{G}_{0}^{0}$ (dashed black line/open symbols) obtained by subspace identification: frequency response (a) gain and (b) phase (the shaded region corresponds to $|\text{arg}(G)-\text{arg}(G_{0}^{0})|\leqslant \unicode[STIX]{x03C0}$ ); (c) poles $\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FF}$ .

We note in passing the monotonically decreasing phase (except for ‘upward jumps’ corresponding to unstable modes) with a nearly constant negative slope in (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 $\unicode[STIX]{x1D705}U$ in the shear layer. Fitting the data with a straight line in the range $5\leqslant \unicode[STIX]{x1D714}\leqslant 25$ leads to a value of $\unicode[STIX]{x1D705}\approx 0.73$ , consistent with the mean flow in figure 6(a). This value of $\unicode[STIX]{x1D705}$ is larger than the value 0.53 obtained by Barbagallo et al. (Reference Barbagallo, Sipp and Schmid2009) 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 $\boldsymbol{B}$ is localized within the incoming boundary layer at $x_{2}^{0}=0.02>0$ and not directly at the wall.

4.2 Structured ${\mathcal{H}}_{\infty }$ -synthesis

After designing a ROM, we seek a controller $K_{1}^{\prime }$ outputting the forcing signal

(4.2) $$\begin{eqnarray}u_{1}^{\prime }:=K_{1}^{\prime }y\end{eqnarray}$$

from the measurement $y$ . When the plant $\widetilde{G}_{0}^{0}$ is interconnected with $K_{1}^{\prime }$ , the closed loop is characterized by four transfer functions

(4.3) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}y\\ u_{1}^{\prime }\end{array}\right]=\left[\begin{array}{@{}cc@{}}\widetilde{G}_{0}^{0}S & S\\ T & K_{1}^{\prime }S\end{array}\right]\left[\begin{array}{@{}c@{}}n\\ w\end{array}\right],\end{eqnarray}$$

linking the two inputs $n$ and $w$ to the two outputs $y$ and $u_{1}^{\prime }$ . These four transfers are combinations of only two interconnected blocks $\widetilde{G}_{0}^{0}$ and $K_{1}^{\prime }$ , therefore they are not independent of each other. In particular, the sensitivity function $S$ and complementary sensitivity function $T$ , respectively defined as

(4.4a,b ) $$\begin{eqnarray}S={\displaystyle \frac{1}{1-\widetilde{G}_{0}^{0}K_{1}^{\prime }}}\quad \text{and}\quad T={\displaystyle \frac{\widetilde{G}_{0}^{0}K_{1}^{\prime }}{1-\widetilde{G}_{0}^{0}K_{1}^{\prime }}},\end{eqnarray}$$

satisfy the relation $S-T=1$ . The controller $K_{1}^{\prime }$ 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^{\prime }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 $\unicode[STIX]{x1D70C}$ bounding the ${\mathcal{H}}_{\infty }$ -norm of the four weighted transfer functions

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \Vert W_{S}S\Vert _{\infty }\leqslant \unicode[STIX]{x1D70C}, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \Vert W_{T}T\Vert _{\infty }\leqslant \unicode[STIX]{x1D70C}, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle \Vert W_{GS}\widetilde{G}_{0}^{0}S\Vert _{\infty }\leqslant \unicode[STIX]{x1D70C}, & \displaystyle\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle \Vert W_{K^{\prime }S}K_{1}^{\prime }S\Vert _{\infty }\leqslant \unicode[STIX]{x1D70C}, & \displaystyle\end{eqnarray}$$

while simultaneously ensuring that the closed-loop transfers and the controller $K_{1}^{\prime }$ remain stable. We recall that for a single-input single-output (SISO) system, the ${\mathcal{H}}_{\infty }$ -norm of a stable transfer function $H$ is given by the maximum of the modulus of the frequency response

(4.9) $$\begin{eqnarray}\Vert H\Vert _{\infty }=\max _{\unicode[STIX]{x1D714}\in \mathbb{R}}|H(\text{i}\unicode[STIX]{x1D714})|,\end{eqnarray}$$

hence, conditions (4.5)–(4.8) are equivalent to

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle |S(\text{i}\unicode[STIX]{x1D714})|\leqslant \unicode[STIX]{x1D70C}/|W_{S}(\text{i}\unicode[STIX]{x1D714})|, & \displaystyle\end{eqnarray}$$

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle |T(\text{i}\unicode[STIX]{x1D714})|\leqslant \unicode[STIX]{x1D70C}/|W_{T}(\text{i}\unicode[STIX]{x1D714})|, & \displaystyle\end{eqnarray}$$

(4.12) $$\begin{eqnarray}\displaystyle & \displaystyle |\widetilde{G}_{0}^{0}S(\text{i}\unicode[STIX]{x1D714})|\leqslant \unicode[STIX]{x1D70C}/|W_{GS}(\text{i}\unicode[STIX]{x1D714})|, & \displaystyle\end{eqnarray}$$

(4.13) $$\begin{eqnarray}\displaystyle & \displaystyle |K_{1}^{\prime }S(\text{i}\unicode[STIX]{x1D714})|\leqslant \unicode[STIX]{x1D70C}/|W_{KS}(\text{i}\unicode[STIX]{x1D714})|, & \displaystyle\end{eqnarray}$$

for all real values of $\unicode[STIX]{x1D714}$ . The design of controller $K_{1}^{\prime }$ 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 $\unicode[STIX]{x1D70C}$ is minimized. If one of the closed-loop constraints is not met, the free parameters are adjusted and $\unicode[STIX]{x1D70C}$ optimized again, in an iterative fashion until all constraints are satisfied.

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}^{\prime }S$ (Henning & King Reference Henning and King2007; Henning et al. Reference Henning, Pastoor, King, Noack, Tadmor and King2007; Williams et al. Reference Williams, Kerstens, Pfeiffer, King and Colonius2010), but not $\widetilde{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 $\widetilde{G}_{0}^{0}K_{1}^{\prime }$ , as done in the loop-shaping approach (Dahan et al. Reference Dahan, Morgans and Lardeau2012; Jones et al. Reference Jones, Heins, Kerrigan, Morrison and Sharma2015; Li & Morgans Reference Li and Morgans2016; Dalla Longa et al. Reference Dalla Longa, Morgans and Dahan2017). Loop shaping takes advantage of the linearity of the loop gain with respect to the controller $K_{1}^{\prime }$ 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 ${\mathcal{H}}_{\infty }$ 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}^{\prime }$ 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 ${\mathcal{H}}_{\infty }$ -synthesis problem is known to be NP-hard (Blondel & Tsitsiklis Reference Blondel and Tsitsiklis1997) 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 Apkarian & Noll (Reference Apkarian and Noll2006).

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 $\unicode[STIX]{x1D70C}$ to a final value of 0.6.

Figure 12.

Frequency response of the four closed-loop transfer functions of the ROM (black lines): (a) $\widetilde{G}_{0}^{0}S$ , (b) $S$ , (c) $T$ , (d) $K_{1}^{\prime }S$ . The frequency templates $\unicode[STIX]{x1D70C}/|W_{GS}|$ , $\unicode[STIX]{x1D70C}/|W_{S}|$ , $\unicode[STIX]{x1D70C}/|W_{T}|$ are plotted in blue ( $\unicode[STIX]{x1D70C}/|W_{K^{\prime }S}|$ is not shown in (d) as it is just a constant curve located well above the maximum of $|K_{1}^{\prime }S|$ ), with solid lines for $\unicode[STIX]{x1D70C}=1$ and dashed lines for the optimized value of the multiplicative constant $\unicode[STIX]{x1D70C}=0.6$ . In (b), the shaded region is the forbidden area where the modulus margin constraint (4.20) is violated.

4.2.1 Shaping $\widetilde{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 $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{h}$ of $G_{0}^{0}(\text{i}\unicode[STIX]{x1D714})$ , we must ensure that

(4.14) $$\begin{eqnarray}|S(\text{i}\unicode[STIX]{x1D714}_{h})|<1,\end{eqnarray}$$

such that the closed-loop gain

(4.15) $$\begin{eqnarray}|\widetilde{G}_{0}^{1}|=|\widetilde{G}_{0}^{0}S|\end{eqnarray}$$

be smaller than the open-loop gain $|\widetilde{G}_{0}^{0}|$ at this frequency. This is achieved by designing a ‘bell-shaped’ frequency template

(4.16) $$\begin{eqnarray}1/W_{GS}=1/k_{GS}{\displaystyle \frac{s^{2}/\unicode[STIX]{x1D6FC}^{2}+2\unicode[STIX]{x1D709}\unicode[STIX]{x1D714}_{h}s+\unicode[STIX]{x1D714}_{h}^{2}/\unicode[STIX]{x1D6FC}^{2}}{s^{2}+2\unicode[STIX]{x1D709}\unicode[STIX]{x1D714}_{h}s+\unicode[STIX]{x1D714}_{h}^{2}}},\end{eqnarray}$$

with $\unicode[STIX]{x1D6FC}^{2}>1$ (see figure 12 a) applied to the closed-loop transfer function $\widetilde{G}_{0}^{0}S$ . The gain $|1/W_{GS}(\text{i}\unicode[STIX]{x1D714})|$ varies between a minimum value of $1/(\unicode[STIX]{x1D6FC}^{2}k_{GS})$ , reached asymptotically when $\unicode[STIX]{x1D714}\rightarrow 0$ and $\unicode[STIX]{x1D714}\rightarrow \infty$ , and a maximum value of $1/k_{GS}$ , reached at the resonance frequency $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{h}$ . The parameter $\unicode[STIX]{x1D709}$ determines the width of the filter. By choosing

(4.17) $$\begin{eqnarray}k_{GS}=1/|\widetilde{G}_{0}^{0}(\text{i}\unicode[STIX]{x1D714}_{h})|\end{eqnarray}$$

and constraining $\widetilde{G}_{0}^{0}S$ , we indirectly enforce a constraint on the sensitivity function

(4.18) $$\begin{eqnarray}|S(\text{i}\unicode[STIX]{x1D714}_{h})|<\unicode[STIX]{x1D70C}.\end{eqnarray}$$

Since the optimal value of $\unicode[STIX]{x1D70C}$ is 0.6 in our case, we enforce a reduction of $60\,\%$ in the height of the main resonance peak at $\unicode[STIX]{x1D714}_{h}$ .

Shaping $\widetilde{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 $\unicode[STIX]{x1D6FC}$ with a zero of $K_{1}^{\prime }$ is not a robust way to suppress the resonance at $\unicode[STIX]{x1D714}_{h}$ since the mean-flow model is only an approximation of the true plant. Worse still, the unstable pole $\unicode[STIX]{x1D6FC}$ will remain present in the closed loop through the transfer $\widetilde{G}_{0}^{0}S$ . A much better alternative is to strongly damp the growth rate associated with this pole by shaping $\widetilde{G}_{0}^{0}S$ instead of $S$ . When targeting the sensitivity function, as in the mixed-sensitivity approach, pole cancellation would typically occur since $S$ is a function of the product $\widetilde{G}_{0}^{0}K_{1}^{\prime }$ , i.e. the loop gain, only. The same problem occurs in the loop-shaping framework, which directly targets the loop gain. When targeting $\widetilde{G}_{0}^{0}S$ instead, pole cancellation does not occur since that transfer depends on both the loop gain and $\widetilde{G}_{0}^{0}$ . The four-block approach is therefore well suited to the control of unstable plants.

In fact, with appropriate choices for $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D709}$ , we see in figure 12(a) that we manage to suppress all three resonance peaks at $\unicode[STIX]{x1D714}_{h}$ and $\unicode[STIX]{x1D714}_{h}\pm \unicode[STIX]{x1D714}_{l}$ , associated with modes $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ . As a result, we verify in (b) that the gain of the sensitivity function is less than 1 at these three frequencies. We also verify in this panel that $\unicode[STIX]{x1D70C}=|S(\text{i}\unicode[STIX]{x1D714}_{h})|=0.6$ . We will later verify in figure 17 from § 5.3 that the corresponding poles $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ 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-\widetilde{G}_{0}^{0}K_{1}^{\prime }$ . 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)

(4.19) $$\begin{eqnarray}\unicode[STIX]{x0394}M=\min _{\unicode[STIX]{x1D714}\in \mathbb{R}}|1-\widetilde{G}_{0}^{0}(\text{i}\unicode[STIX]{x1D714})K_{1}^{\prime }(\text{i}\unicode[STIX]{x1D714})|\end{eqnarray}$$

be sufficiently large. In practice, we ask that $\unicode[STIX]{x0394}M>-6$  dB, or $\unicode[STIX]{x0394}M>0.5$ in linear scale, which amounts to bounding the sensitivity function

(4.20) $$\begin{eqnarray}\Vert S\Vert _{\infty }<2.\end{eqnarray}$$

The constraint on $S(\text{i}\unicode[STIX]{x1D714})$ being independent of $\unicode[STIX]{x1D714}$ , we consider a constant frequency template $1/W_{S}$ and optimize $\unicode[STIX]{x1D70C}$ such that

(4.21) $$\begin{eqnarray}|S(\text{i}\unicode[STIX]{x1D714})|\leqslant {\displaystyle \frac{\unicode[STIX]{x1D70C}}{|W_{S}|}}<2\end{eqnarray}$$

for any real $\unicode[STIX]{x1D714}$ . In figure 12(b), constraint (4.20) is represented by the shaded region. We see that the value of $\unicode[STIX]{x1D70C}=0.6$ reached after optimization did not allow the constraint to be satisfied in the vicinity of $\unicode[STIX]{x1D714}\approx 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}$ .

4.2.3 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}^{\prime }(0)=0$ , is equivalent to a condition on the complementary sensitivity function

(4.22) $$\begin{eqnarray}T(0)=0,\end{eqnarray}$$

which can be approximately achieved using constraint (4.11) and a filter $W_{T}$ with a sufficiently large static gain $|W_{T}(0)|\gg 1$ , such that

(4.23) $$\begin{eqnarray}|T(0)|<{\displaystyle \frac{\unicode[STIX]{x1D70C}}{|W_{T}(0)|}}\ll 1.\end{eqnarray}$$

The frequency template $1/W_{T}$ is therefore defined as a high-pass filter, as can be verified in figure 12(c).

5.1 Time scales

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}^{\prime }$ is added to the loop, its internal state is initialized to 0. The internal states of the other controllers $K_{j<m}^{\prime }$ 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}^{\prime }$ . 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 time scale of the order of $O(1)$ , which seems sufficient for the system to latch onto a new state. This time scale 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 time scales in the transient evolution is once again related to the two physical processes occurring 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, time scale $O(0.1)$ , highlighted in very dark grey and visible at iteration 1 only. This time scale is intrinsic to the controller and the large ‘overshoot’ observed at the first iteration is due to the transition from a 2-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}$ .

5.2 Attractors

In figure 8 from § 3, we plotted a spectrogram of the time series $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 2-torus characterized by the two fundamental frequencies $\unicode[STIX]{x1D714}_{l}$ and $\unicode[STIX]{x1D714}_{h}$ . But soon after $t_{1}$ , the flow transitions to a limit cycle with a fundamental frequency close to the previous harmonic $\unicode[STIX]{x1D714}_{h}+\unicode[STIX]{x1D714}_{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 2-torus, which is strongly dominated by its fast fundamental frequency. The trace of interactions between the high-frequency fundamental and a non-zero low-frequency mode at $\unicode[STIX]{x1D714}\approx 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 (a). The decrease of the oscillation amplitude at the sensor $y$ is less dramatic (b). In figure 14, we project each attractor in the phase plane $[\unicode[STIX]{x0394}u_{2}(1/4,0),\unicode[STIX]{x0394}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.

Figure 14.

Characterization of dynamical equilibria $m=0,\ldots ,5$ in the phase plane $[\unicode[STIX]{x0394}u_{2}(1/4,0),\unicode[STIX]{x0394}u_{2}(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 2-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 2-torus: the extra dimension of the attractor is so thin that it is hidden by the line width.

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.

Figure 15.

Vorticity snapshots at a time where the shear-stress measurement $y$ is maximum (a) $m=0$ , (b) $m=1$ , (c) $m=2$ , (d) $m=3$ , (e) $m=4$ , (f)  $m=5$ . The colour map has been cropped to the same arbitrary extremal values in all plots for comparison.

Figure 16.

Evolution of the gain $|G|$ during the iterative procedure. (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-1}^{m}$ (dashed lines) at $t=t_{m}^{+}$ to $G_{m}^{m}$ (solid lines) at $t=t_{m+1}^{-}$ (with $t_{6}\rightarrow \infty$ ), under the effects of nonlinearity. Same labels as figure 15. Vertical lines indicate Koopman frequencies of the permanent state obtained with the controller $K_{m}$ (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.

5.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 developed regime. The plant $G_{m}^{m}$ (solid lines), computed using (2.32), characterizes the fully developed flow in feedback loop with the controller $K_{m}$ . At $t=t_{m+1}$ , a new controller correction $K_{m+1}^{\prime }$ 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. (4.16)) centred about the dominant Koopman mode. This procedure leads to the shaped plant $G_{m}^{m+1}$ (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_{5}^{4}$ is therefore nearly equal to the last 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}^{m+1}$ are indicated with open symbols. In fact, the poles are computed from a ROM $\widetilde{G}$ obtained by subspace identification of $G$ , at each step. We only follow the evolution of the four poles $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ , $\unicode[STIX]{x1D6FE}$ , $\unicode[STIX]{x1D6FF}$ 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}\rightarrow G_{m}^{m+1}$ , the plant is stabilized. The growth rate of all four resonant poles is damped (except pole $\unicode[STIX]{x1D6FE}$ 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 with the frequency of the dominant Koopman mode, specifically targeted by the ‘bell-shaped’ template. Then, during the relaxation phase $G_{m}^{m+1}\rightarrow 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 developed 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 $\{\text{flow}+\text{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 $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ , $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FF}$ , which synchronize with mode $b$ . Then the controlled flow switches to a limit cycle driven by mode $\unicode[STIX]{x1D6FD}$ , 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 $\unicode[STIX]{x1D6FC}$ , whereas it is associated with mode $\unicode[STIX]{x1D6FD}$ when $K_{1}^{\prime }$ is in the loop, which explains the sudden frequency jump in the spectrogram. Interestingly, we note the complete disappearance of the low frequency at $\unicode[STIX]{x1D714}_{l}$ , despite the fact that controller $K_{1}^{\prime }$ did not target the corresponding mode $b$ . The suppression of this oscillation probably occurred through a nonlinear mechanism: by damping modes $\unicode[STIX]{x1D6FD}$ , $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FF}$ corresponding to the oscillations at $\unicode[STIX]{x1D714}_{h}-\unicode[STIX]{x1D714}_{l}$ , $\unicode[STIX]{x1D714}_{h}+\unicode[STIX]{x1D714}_{l}$ and $\unicode[STIX]{x1D714}_{h}+2\unicode[STIX]{x1D714}_{l}$ , the controller $K_{1}^{\prime }$ has indirectly targeted the low frequency $\unicode[STIX]{x1D714}_{l}$ . In subsequent iterations $m>1$ , the flow remains on a limit cycle governed by mode $\unicode[STIX]{x1D6FD}$ , which explains the smooth evolution of the dominant frequency in the spectrogram (figure 8), and the qualitatively similar phase portraits in figure 14.