2.1 In-flow control set-up

In the first control set-up, the momentum equation (2.2) is forced by a pair of antisymmetric body forces $\boldsymbol{f}^{\prime }$ that serve as an in-flow actuator:

(2.3)$$\begin{eqnarray}\boldsymbol{f}^{\prime }=\boldsymbol{B}(x,y)[u(t)+w(t)],\end{eqnarray}$$

where $u(t)$ is the actuator signal provided by a transverse velocity sensor positioned a distance $d$ downstream of the cylinder, $w(t)$ is a disturbance from the actuator and $\boldsymbol{B}(x,y)$ is the spatial distribution of the actuator and disturbance. More specifically,

(2.4)$$\begin{eqnarray}\boldsymbol{B}=\left[\begin{array}{@{}c@{}}\cos \unicode[STIX]{x1D703}({\mathcal{S}}(A,\unicode[STIX]{x1D70E},r,\unicode[STIX]{x1D703})-{\mathcal{S}}(A,\unicode[STIX]{x1D70E},r,-\unicode[STIX]{x1D703}))\\ \sin \unicode[STIX]{x1D703}({\mathcal{S}}(A,\unicode[STIX]{x1D70E},r,\unicode[STIX]{x1D703})+{\mathcal{S}}(A,\unicode[STIX]{x1D70E},r,-\unicode[STIX]{x1D703})),\end{array}\right],\end{eqnarray}$$

where

(2.5)$$\begin{eqnarray}{\mathcal{S}}(A,\unicode[STIX]{x1D70E},r,\unicode[STIX]{x1D703})=\frac{A}{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}^{2}}\exp \left(-\frac{(x-r\cos \unicode[STIX]{x1D703})^{2}+(y-r\sin \unicode[STIX]{x1D703})^{2}}{2\unicode[STIX]{x1D70E}^{2}}\right),\end{eqnarray}$$

where the magnitude is $A=1.0$ and the standard deviation $\unicode[STIX]{x1D70E}=0.1$. The centres of the above distributions are near the separation points at a distance of $r=0.6$ from the cylinder’s centre and an angle of $\unicode[STIX]{x1D703}=\pm 70^{\circ }$ from the cylinder’s downstream-pointing horizontal. We fix the form and the location of the actuator and vary only the sensor position along the centreline to investigate the influence of sensor placement. This arrangement is similar to that used by Illingworth (Reference Illingworth2016).

2.2 Body-mounted control set-up

The second control set-up considered uses a body-mounted actuator and a body-mounted sensor. Similar to a fluid–structure interaction system, we consider an oscillatory cylinder whose transverse acceleration $a$ is controlled according to the feedback signal provided by a lift sensor attached on the cylinder. Optimal feedback controllers are designed for different Reynolds numbers to suppress vortex shedding and the vibration of the bluff body.

Figure 1.

(a) A pair of body forces $\boldsymbol{f}^{\prime }$ applied near the cylinder surface cooperates with a velocity sensor placed on the centreline at a distance $d$ from the centre of the cylinder. (b) The force sensor is placed on the cylinder to measure the lift $l^{\prime }$ which is fed to the actuator that controls the transverse acceleration $a$ of the cylinder.

3.1 State-space formulation and system identification

The linearised Navier–Stokes equations can be written in standard state-space form which is useful for performing input–output analysis. To evaluate the linear dynamics of a flow system, it is convenient to take Laplace transforms. Introducing the transformation into (2.2), we obtain the equations

(3.1a,b)$$\begin{eqnarray}s\hat{\boldsymbol{u}}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\hat{\boldsymbol{u}}+\hat{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U}=-\unicode[STIX]{x1D735}\hat{p}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\hat{\boldsymbol{u}}+\hat{\boldsymbol{f}},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}}=0,\end{eqnarray}$$

where $(\hat{\boldsymbol{u}},~\hat{p})$ and $\hat{\boldsymbol{f}}$ represent the (complex-valued) spatial structure of the velocity and body forcing and $s=\unicode[STIX]{x1D70E}+j\unicode[STIX]{x1D714}$ is the Laplace variable. Thus, the transfer function between the external forcing and response can be written as

(3.2)$$\begin{eqnarray}\left[\begin{array}{@{}c@{}}\hat{\boldsymbol{u}}\\ \hat{p}\end{array}\right]=(s{\mathcal{E}}-{\mathcal{A}})^{-1}\left[\begin{array}{@{}c@{}}\hat{\boldsymbol{f}}\\ 0\end{array}\right],\end{eqnarray}$$

where $(s{\mathcal{E}}-{\mathcal{A}})^{-1}$ is known as the resolvent operator and ${\mathcal{A}}$ is the linearised Navier–Stokes operator around the base flow:

(3.3a,b)$$\begin{eqnarray}{\mathcal{A}}=\left[\begin{array}{@{}cc@{}}-\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}-()\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2} & -\unicode[STIX]{x1D735}\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }() & 0\end{array}\right],\quad {\mathcal{E}}=\left[\begin{array}{@{}cc@{}}I & 0\\ 0 & 0\end{array}\right].\end{eqnarray}$$

We now consider writing (3.2) in state-space form for a linear, time-variant dynamical system $P(s)$:

$$\begin{eqnarray}s{\mathcal{E}}\boldsymbol{x}={\mathcal{A}}\,\boldsymbol{x}+{\mathcal{B}}e,\end{eqnarray}$$

(3.4)$$\begin{eqnarray}y={\mathcal{C}}\boldsymbol{x}+{\mathcal{D}}e,\end{eqnarray}$$

where $\boldsymbol{x}=[\hat{\boldsymbol{u}}~~\hat{p}]^{\text{T}}$ is the system state, $e$ is an input vector of dimension $p$ and $y$ is an output vector of dimension $q$. The feedback control arrangement used here is single-input single-output with $p=q=1$. The vector ${\mathcal{B}}$ is determined by the shape of the actuation, which is from the spatial discretisation of the external forcing term $\hat{\boldsymbol{f}}$. The matrices ${\mathcal{C}}$ and ${\mathcal{D}}$ represent output and feed-forward dynamics, respectively. The transfer function between the output $y$ and input $e$ can then be written as

(3.5)$$\begin{eqnarray}P(s)={\mathcal{C}}(s{\mathcal{E}}-{\mathcal{A}})^{-1}{\mathcal{B}}+{\mathcal{D}}.\end{eqnarray}$$

In general, $P(s)$ is of high dimension, which makes the control design problem computationally intractable. However, it is feasible to solve this linear system and get the frequency response data from the actuation to the measurement. Reduced-order models can then be formed for the corresponding frequency responses using the vector-fitting algorithm VECTFIT (Gustavsen & Semlyen Reference Gustavsen and Semlyen1999; Gustavsen Reference Gustavsen2006; Deschrijver et al. Reference Deschrijver, Mrozowski, Dhaene and De Zutter2008). VECTFIT identifies a ROM $\widetilde{P}(s)$ of a significantly smaller dimension $N$ in the rational form

(3.6)$$\begin{eqnarray}\widetilde{P}(s)=\mathop{\sum }_{m=1}^{N}{\displaystyle \frac{r_{m}}{s-a_{m}}}+sh+d,\end{eqnarray}$$

where the poles $a_{m}$, residues $r_{m}$ and terms $d$ and $h$ are identified such that the frequency response of $\widetilde{P}(s)$ is within a distance $\unicode[STIX]{x1D716}$ of the original system $P(s)$ over a broad frequency range.

3.2 Feedback controller design

The feedback controller is designed based on the ROM $\widetilde{P}(s)$ using ${\mathcal{H}}_{\infty }$ loop shaping. The control block diagram is shown in figure 2, where ${\hat{w}}$ represents disturbances from the actuator and $n$ represents noise at the sensor. We assume a positive feedback configuration as illustrated in figure 2. The controller $K(s)$ for plant $P(s)$ is designed using the loop-shaping design method of Glover & McFarlane (Reference Glover and McFarlane1989) which maximises the normalised coprime stability margin $b(\widetilde{P},K)$ of the plant–controller feedback loop

(3.7)$$\begin{eqnarray}b=\left\Vert \left[\begin{array}{@{}c@{}}K\\ I\end{array}\right](I+\widetilde{P}K)^{-1}\left[\begin{array}{@{}cc@{}}\widetilde{P} & I\end{array}\right]\right\Vert _{\infty }^{-1},\end{eqnarray}$$

where $b\in [0,1]$, which can be maximised over all stabilising $K$ to give

(3.8)$$\begin{eqnarray}b_{opt}(\widetilde{P})=\sup _{K}b(\widetilde{P},K).\end{eqnarray}$$

Physically, $b_{opt}$ is an indication of the robustness of the closed-loop system to unmodelled dynamics and also serves as a performance measure: the larger the value of $b_{opt}$, the better the performance and the greater the robustness of the closed-loop system.

For single-input single-output systems as considered here, a compensator which weights the plant according to the control objectives is used of the form

(3.9)$$\begin{eqnarray}W(s)=k\frac{a^{2}}{(s+a)^{2}}.\end{eqnarray}$$

The parameters $k$ and $a$ are chosen such that the gain of the weighted plant is sufficiently high at frequencies where good disturbance attenuation is required and is sufficiently low at high frequencies where modelling uncertainties will be greatest (Illingworth Reference Illingworth2016).

Figure 2.

Block diagram of the closed-loop control system.

3.3 Numerical set-up

We consider an incompressible flow with free-stream velocity $U_{\infty }$ past a two-dimensional circular cylinder of diameter $D$. Simulations are conducted on the computing platform FEniCS (Logg, Mardal & Wells Reference Logg, Mardal and Wells2012) which has been extensively used for fluid mechanics (Mortensen, Langtangen & Wells Reference Mortensen, Langtangen and Wells2011; Nguyen et al. Reference Nguyen, Jansson, Goude and Hoffman2019; Vasilyeva et al. Reference Vasilyeva, Chung, Efendiev and Kim2019). Three versions of the flow are solved: the steady solution of (2.1) (base flow), the linearised perturbation equation (3.4) in state-space form and the fully nonlinear Navier–Stokes equations. We employ the same computational domain as that used by Leontini et al. (Reference Leontini, Stewart, Thompson and Hourigan2006), as shown in figure 3, and discretise using Taylor–Hood finite elements over a structured mesh. In order to appropriately resolve the details of the flow, the mesh points are clustered smoothly near the cylinder and in the wake. More specifically, the mesh consists of $5.46\times 10^{4}$ triangles and the minimum wall-normal size around the cylinder is 0.01. A backward Euler scheme is used for time discretisation ($\unicode[STIX]{x0394}t=0.01$) in the direct numerical simulations with the maximum Courant number below 0.6 to ensure accuracy.

Figure 3.

Computational domain and boundary conditions for the steady and unsteady linearised Navier–Stokes equations for the in-flow set-up.

The base flow $(\boldsymbol{U},P)$, which is governed by the unforced steady Navier–Stokes equations, is solved using a Newton method. The corresponding boundary conditions are summarised in figure 3. Dirichlet boundary conditions are imposed at the inlet $\unicode[STIX]{x1D6E4}_{in}~(x<0)$ and at the cylinder surface $\unicode[STIX]{x1D6E4}_{wall}$ (centred at $x=y=0$):

(3.10a,b)$$\begin{eqnarray}\displaystyle U=1,\quad V=0\quad \text{on }\unicode[STIX]{x1D6E4}_{in}, & & \displaystyle\end{eqnarray}$$

(3.10c,d)$$\begin{eqnarray}\displaystyle U=0,\quad V=0\quad \text{on }\unicode[STIX]{x1D6E4}_{wall}. & & \displaystyle\end{eqnarray}$$

Symmetry conditions are enforced on the top and bottom boundaries $\unicode[STIX]{x1D6E4}_{top}$ and $\unicode[STIX]{x1D6E4}_{bottom}$ ($0<x/D<23,y/D=\pm 15$):

(3.11a,b)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}}=0,\quad V=0\quad \text{on }\unicode[STIX]{x1D6E4}_{top}\cup \unicode[STIX]{x1D6E4}_{bottom},\end{eqnarray}$$

while the pressure and the velocity are combined into the standard outflow conditions on the outlet boundary $\unicode[STIX]{x1D6E4}_{out}$ ($x/D=23,-15<y/D<15$):

(3.12)$$\begin{eqnarray}-P\boldsymbol{n}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D735}\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{n}=0\quad \text{on }\unicode[STIX]{x1D6E4}_{out},\end{eqnarray}$$

where $\boldsymbol{n}$ denotes the outward-pointing normal vector on the boundary. The base flow is the same for all control set-ups. Vorticity contours for the base flow are shown in figure 4(a) at $Re=60,80$ and $100$.

The basic boundary conditions for the perturbation $(\boldsymbol{u}^{\prime },p^{\prime })$ are also depicted in figure 3. For the in-flow control set-up, the linear perturbation system has the same boundary conditions as the base flow except at the inlet where homogeneous boundary conditions $(u^{\prime }=v^{\prime }=0)$ are enforced to ensure zero perturbations at infinity. In the body-mounted control set-up, the flow is solved in an accelerated frame of reference attached to the cylinder instead of moving the cylinder directly. For this, the transverse acceleration of the frame $a$ is treated as an extra forcing term $\boldsymbol{f}^{\prime }=[0~~a]^{\text{T}}$. Thus, the boundary conditions at the inlet, top and bottom boundaries are modified (Leontini et al. Reference Leontini, Stewart, Thompson and Hourigan2006):

(3.13a,b)$$\begin{eqnarray}\displaystyle u^{\prime }=0,\quad v^{\prime }=\int _{0}^{t}a(t)\,\text{d}t\quad \text{on }\unicode[STIX]{x1D6E4}_{in}, & & \displaystyle\end{eqnarray}$$

(3.13c,d)$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}y}}=0,\quad v^{\prime }=\int _{0}^{t}a(t)\,\text{d}t\quad \text{on }\unicode[STIX]{x1D6E4}_{top}\cup \unicode[STIX]{x1D6E4}_{bottom}, & & \displaystyle\end{eqnarray}$$

where zero initial conditions are assumed. In the linear state-space model (3.4), these boundary conditions are enforced in the frequency domain by applying Laplace transforms:

(3.14a,b)$$\begin{eqnarray}\displaystyle \hat{u} =0,\quad \hat{v}={\displaystyle \frac{\hat{a}}{s}}\quad \text{on }\unicode[STIX]{x1D6E4}_{in}, & & \displaystyle\end{eqnarray}$$

(3.14c,d)$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\hat{u} }{\unicode[STIX]{x2202}y}}=0,\quad \hat{v}={\displaystyle \frac{\hat{a}}{s}}\quad \text{on }\unicode[STIX]{x1D6E4}_{top}\cup \unicode[STIX]{x1D6E4}_{bottom}. & & \displaystyle\end{eqnarray}$$

The stability analysis of the base flows and discretised perturbation systems have been validated using the results of Barkley (Reference Barkley2006). To validate the controllers, direct numerical simulations are performed using the IPCS (Incremental Pressure Correction Scheme) method which has been extensively tested in Logg et al. (Reference Logg, Mardal and Wells2012). The corresponding boundary conditions are modified based on the identity $(\boldsymbol{u},p)=(\boldsymbol{U}+\boldsymbol{u}^{\prime },~P+p^{\prime })$.

4.1 Open-loop system identification

Figure 4(b) shows the normalised vorticity field of the corresponding perturbation system $P(s)$ actuated by the body force $\hat{\boldsymbol{f}}$ at the unstable frequency $s=j\unicode[STIX]{x1D714}_{u}$ (i.e. the resolvent operator between the input and output in (3.2)). The shear layer emanating from the separation point is disturbed by the actuator and grows into a large-scale vortical structure downstream. As the Reynolds number increases, the large vortical structure formed by the shear layer appears increasingly far downstream. This is because the length of the absolutely unstable region behind the circular cylinder increases with increasing Reynolds number. It is noteworthy that the absolutely unstable region has the opposite trend of length variation for the linear flow (base flow) when compared to the nonlinear (mean) flow (Pier Reference Pier2002).

After obtaining the base flow, the transfer function between the actuator signal and the sensor measurement, denoted by $\widetilde{P}(s)$, is to be determined for control design. Instead of harvesting data from direct numerical simulations (e.g. Dahan, Morgans & Lardeau Reference Dahan, Morgans and Lardeau2012; Illingworth Reference Illingworth2016), we obtain frequency responses by solving (3.5) directly for harmonic inputs ($s=j\unicode[STIX]{x1D714}$) over a broad range of frequencies. The ROMs are identified from the gain and phase of the response data by utilising a vector-fitting algorithm (Gustavsen Reference Gustavsen2013). The frequency range chosen here is $\unicode[STIX]{x1D714}\in [10^{-6},10^{6}]$  rad s$^{-1}$ with 250 sampling points spaced logarithmically. At each Reynolds number, the computational cost is less than five minutes using a direct LU solver (Amestoy et al. Reference Amestoy, Duff, Koster and L’Excellent2001, Reference Amestoy, Guermouche, L’Excellent and Pralet2006) and a standard node on a cluster. Figure 5(a) shows comparisons of response data near the unstable frequency ($\unicode[STIX]{x1D714}\approx 0.8~\text{rad}~\text{s}^{-1}$) between the linear perturbation systems $P(s)$ and the identified ROMs $\widetilde{P}(s)$ at $Re=60,80$ and $100$ with the velocity sensor placed at $d=2.8D$. The order (or dimension) of $\widetilde{P}(s)$ is chosen such that the fitting residual $\unicode[STIX]{x1D716}$ is below 10^-5. This is achieved with an order of 30 or less for all Reynolds numbers considered.

The accuracy of the ROMs is validated by the excellent agreement between the open-loop impulse responses from these models and direct numerical simulations, which is summarised in figure 5(b). Due to the small-perturbation assumption, the magnitude of the impulse equals $1\times 10^{-4}$ and the linearity of the flow systems is also confirmed. Thus, the ROMs found are good approximations of the true systems.

Figure 5.

(a) Frequency responses from $P(j\unicode[STIX]{x1D714})$ (——) compared to those from ROMs $\widetilde{P}(j\unicode[STIX]{x1D714})$ ($\bullet$, blue) at $Re=60,80$ and $100$. (b) The corresponding open-loop impulse responses from numerical simulations. The results for $Re=60$ and $80$ are multiplied by 20 and 3, respectively, so that the same scale can be used.

4.2 Model-based feedback control

Based on the identified ROMs, we design an optimal feedback controller for each model and ensure that the controller satisfies closed-loop stability. The parameters in the compensator weight are determined after searching different value combinations to achieve the largest optimal stability margin $b_{opt}$. The controllers are implemented in the full nonlinear Navier–Stokes system with an impulse from the actuator of magnitude 10^-3. The resulting performance of the controllers for each of the cases described in § 4.1 is shown in figure 6.

Figure 6.

Direct numerical simulation results of closed-loop systems. (a) Time evolution of the transverse velocity at the sensor ($\bullet$) and the total perturbation energy $E(t)$ in log scale at $Re=60$ (– ⋅ –, blue), $Re=80$ (- -, red) and $Re=100$ (——). (b) Vorticity contours (dashed lines for negative and solid lines for positive vorticity) for the perturbation systems at $t=75$ (▴) at $Re=60,80$ and $100$ (from top to bottom). All contour plots share the same colour range. (c) Table of parameters.

The stabilising effect of the controller for each case is observed from the time evolution of the transverse velocity ($v_{1}^{\prime }/v_{2}^{\prime }/v_{3}^{\prime }$) at the sensor, as shown in figure 6(a). More convincing evidence for the complete suppression of vortex shedding is the total perturbation energy plotted beneath, which is defined as

(4.1)$$\begin{eqnarray}E(t)={\displaystyle \frac{1}{2}}\int _{\unicode[STIX]{x1D6FA}}({u^{\prime }}^{2}(x,~y,~t)+{v^{\prime }}^{2}(x,~y,~t))\,\text{d}x\,\text{d}y,\end{eqnarray}$$

where $u^{\prime }$ and $v^{\prime }$ are the streamwise and transverse perturbation velocity components, respectively. It is clear that the controller for $Re=60$ performs best with the strongest attenuation. Although the perturbation energy for $Re=100$ remains large after $t=150$, a decreasing trend is quite clear. Figure 6(b) shows the instantaneous vorticity perturbation fields for all three cases at $t=75$. The perturbation vorticity at $Re=60$ stays small without observable vortex shedding whereas at a higher Reynolds number of 80, the vorticity remains small near the sensor but is stronger downstream. At the highest Reynolds number of 100, the vorticity is far stronger with clear vortex shedding, which is consistent with the oscillations of $v_{3}^{\prime }$. From the comparison between all three cases, we can see that perturbations are harder to control at higher Reynolds numbers.

The performance of each controller can be quantified using the optimal stability margin $b_{opt}$, which is summarised in figure 6(c) together with other important parameters. Based on the comparison of $b_{opt}$ between the three cases, we can draw a conclusion consistent with the analysis above: the higher the Reynolds number, the smaller the optimal stability margin $b_{opt}$ achieved by the optimal feedback controller.

Figure 7.

(a) Optimal sensor locations (the ridge - -) and contour plot of optimal stability margin $b_{opt}$ against Reynolds number and sensor location $d$. And the largest $b_{opt}$ (——) that can be achieved at different Reynolds numbers is plotted beneath. (b) Loci of unstable poles ($\times$) and critical zeros ($\bullet$, blue; °, red) of transfer functions $\widetilde{P}(s)$ for two cases. Top: different sensor locations at $Re=80$. Bottom: different Reynolds numbers with a sensor placed at $d=2.5D$ downstream of the cylinder.

4.3 Optimal sensor placements

The preliminary investigation summarised in figure 6 indicates a severe deterioration of control performance with increasing Reynolds number for a fixed velocity sensor. To draw more general conclusions about control performance, we vary the position of the velocity sensor along the centreline at different Reynolds numbers. For each case, the position and the form of the actuator are unchanged, and we identify a new ROM from response data at the corresponding Reynolds number. The performance of the ${\mathcal{H}}_{\infty }$-optimal controller for each case is quantified by the optimal stability margin $b_{opt}$ and summarised as a function of Reynolds number and sensor location $d$ in figure 7(a).

First, we focus on the optimal sensor location where the optimal controller shows the best performance at each Reynolds number. Generally, the ideal position for a sensor should allow not only the measurement of the instability developing downstream but also the timely feedback of information to the actuator. Figure 7(a) first shows a contour map of the optimal stability margin $b_{opt}$ against Reynolds number and sensor location $d$. It is clear that a ridge exists which indicates the optimal sensor location as a function of Reynolds number. At each Reynolds number, the optimal controller for the sensor at the ridge line performs better than those for other sensor locations. However, the optimal stability margin $b_{opt}$ on the ridge decreases sharply with increasing Reynolds number, which is shown in the panel beneath the contour map.

Therefore, a fundamental trade-off can be concluded from figure 7(a): the sensor should be close enough to the cylinder to reduce the time delay due to convection, but it should also be far enough from the cylinder to measure important information (e.g. unstable eigenmodes) developing downstream. The compromise between these two conflicting requirements becomes harder to satisfy with increasing Reynolds number, which leads to the optimal sensor location moving downstream linearly. Similar results have been observed in recent work (Oehler & Illingworth Reference Oehler and Illingworth2018) that considers feedback control of the linearised Ginzburg–Landau system.

Generally, the performance and robustness of an optimal controller, as quantified by the optimal stability margin $b_{opt}$, can be linked to the zeros and poles of the corresponding system. Mathematically, poles and zeros of a system are roots of the denominator and numerator of the corresponding transfer function, which determine whether the system is stable, and how the system performs. More specifically, poles capture the form of each component in the system response, whereas zeros reflect how these components combine together, including the phase and magnitude of each component generated by each pole. To investigate these roots in the perturbation system, we consider two cases: (i) fixing the Reynolds number at $Re=80$ and moving the sensor from $d=1D$ to $d=4D$; (ii) fixing the sensor position at $d=2.5D$ and increasing the Reynolds number from 50 to 110. The root loci of these two cases are computed, where unstable poles and two kinds of critical zeros (which we label I and II) are identified and plotted in figure 7(b).

In the first case (top panel in figure 7b), the sensor is moved away from the cylinder, thus measuring information further downstream. This results in the critical zero I moving from the right-half plane (RHP) into the left-half plane (LHP) and the critical zero II moving from the LHP into the RHP. The existence of RHP zeros is problematic for control design because they limit the maximum bandwidth or the maximum frequency that can be controlled with good performance and robustness, as described by Zhou, Doyle & Glover (Reference Zhou, Doyle and Glover1996) and Hoagg & Bernstein (Reference Hoagg and Bernstein2007). In this case, zero I stays in the LHP for $d>1.6D$ whereas zero II stays in the LHP for $d<2.4D$. The optimal sensor location at the Reynolds number considered is $d\approx 2.0D$, for which all zeros stay in the LHP (as indicated by ° (red) in figure 7). This is also consistent with the fundamental trade-off described above.

Increasing Reynolds number also moves the zeros into the RHP, which can be seen from the root loci plotted beneath. In this case, we fix the sensor location at $d=2.5D$ and increase Reynolds number from 50 to 110 in intervals of 10. It is interesting to note that increasing Reynolds number not only increases the real part of the unstable pole but also moves these critical zeros towards the RHP. Based on the root loci in the two cases, we can conclude that it becomes harder to find a good sensor location (where no RHP zeros occur) at higher Reynolds numbers. This difficulty leads to a degradation in the performance and robustness of the optimal controllers at higher Reynolds numbers, as depicted in figure 7(a). Furthermore, zero I appears near the unstable pole and moves closer to the unstable pole at higher Reynolds numbers (or more upstream sensor locations), whereas zero II remains at higher frequencies. Due to the bandwidth limitation from the RHP zeros, the optimal control design algorithm would be able to compute a better controller if zero I stays in the LHP. Thus, the systems would prefer a sensor placed further downstream to prevent RHP zero I. This preference is shown by the contour map in figure 7(a), where a gentle slope of $b_{opt}$ occurs if the sensor is placed downstream of the optimal location but a rapid drop occurs if the sensor is placed upstream.

Similar maps of system roots are also summarised and analysed in the work of Belson et al. (Reference Belson, Semeraro, Rowley and Henningson2013), where optimal controllers were designed for a linearised two-dimensional Blasius boundary layer controlled by different types and positions of actuators and sensors. A degradation of the controllers’ performance and robustness was observed when RHP zeros occurred. Generally, the physical mechanisms behind the RHP zeros are due to (i) the time delay or (ii) the observability of the structures that are to be controlled. In a flow system, when the sensor is far downstream of the actuator, it measures the effect of the actuator with a time delay due to the convective nature of the flow. That is, the sensor measures flow structures that convected past the actuator at an earlier time. With outdated information, the controller poorly estimates and controls the flow structures near the actuator. This time delay becomes more significant as the sensor moves downstream and results in RHP zeros in the reduced-order transfer function $\widetilde{P}(s)$.

However, with a sensor close to the actuator, the performance and robustness of the optimal controller are still restricted by RHP zeros that occur near the unstable pole. Such RHP zeros cancel the effect of the unstable pole and prevent the sensor from measuring the instability. In other words, the poor performance and robustness of the controller are caused by a lack of observability of the unstable mode instead of excessive time delay.

In this section, we have shown that different sensor locations and Reynolds numbers have similar properties that restrict the performance and robustness of the optimal controllers. At higher Reynolds numbers, even an optimal controller performs poorly for both control set-ups. In other words, the best possible performance that can be achieved is severely restricted. From the perspective of control theory, we observe RHP zeros which limit control performance. Thus, we cannot always find a controller with good performance and robustness, and this can be attributed to the compromise between the observability of the instability and the size of the convective time delay.

5.1 Open-loop system identification

The purpose of both control set-ups is to eliminate the perturbations and drive the system towards the steady solution: the base flow. Figure 8 shows the normalised vorticity field of the corresponding perturbation system $P(s)$ actuated by the oscillation of the cylinder at the instability frequency (i.e. the resolvent operator between the input and output in (3.2)). Similar to the in-flow set-up, the large vortical structure actuated on by the moving cylinder develops further downstream as the Reynolds number increases.

The system identification procedure is carried out in a similar manner to § 4.1 and is summarised in figure 9. The ROMs $\widetilde{P}(s)$ are chosen such that the fitting residual $\unicode[STIX]{x1D716}$ is below 10^-5 with orders less than 35. The Bode plots of identified transfer functions between the actuator and the sensor are shown in figure 9(a) and compared to the frequency responses of the true systems. Unlike the in-flow set-up, the perturbation system with body-mounted set-up has infinite zero-frequency response (i.e. the system contains an integrator) and constant infinite-frequency response (i.e. the system contains a non-zero feed-forward term).

Figure 8.

Vorticity contours (dashed lines for negative and solid lines for positive vorticity) for the perturbation system (normalised, real part) actuated by the harmonic oscillation of the cylinder at the unstable frequency for (a) $Re=60$, (b) $Re=80$ and (c) $Re=100$. Contour plots share the same scale.

Figure 9.

(a) Frequency responses from $P(j\unicode[STIX]{x1D714})$ (——) compared to those from ROMs $\widetilde{P}(j\unicode[STIX]{x1D714})$ ($\bullet$, blue) at $Re=60,80$ and $100$. (b) The corresponding open-loop impulse responses from numerical simulations. The results for $Re=60$ and $80$ are multiplied by 15 and 3, respectively, so that the same scale can be used.

Figure 9(b) shows comparisons of open-loop impulse responses (of magnitude 10^-4) from the identified models and direct numerical simulations. The excellent agreement observed validates the accuracy of the ROMs.

5.2 Model-based feedback control

Following the same procedure as in § 4.2, we design optimal controllers for ROMs and implement them in the full nonlinear Navier–Stokes system actuated by an initial impulse of magnitude 10^-4. The parameters of controllers and the corresponding closed-loop simulations are summarised in figure 10. The stabilisation of vortex shedding is achieved only up to $Re=100$, which can be seen both in the time evolution of the lift $l^{\prime }(t)$ and in the total perturbation energy $E(t)$. The comparison among simulations at three Reynolds numbers indicates a similar deterioration of control performance to that seen for the in-flow set-up of § 4.

Figure 10.

Direct numerical simulation results of closed-loop systems. (a) Time evolution of the cylinder lift and the total perturbation energy $E(t)$ in log scale at $Re=60$ (– ⋅ –, blue), $Re=80$ (- -, red) and $Re=100$ (——). (b) Vorticity contours (dashed lines for negative and solid lines for positive vorticity) for the perturbation systems at $t=75$ (▴) at $Re=60,80$ and $100$ (from top to bottom). All contour plots share the same colour range. (c) Table of parameters.

Figure 11.

(a) The largest $b_{opt}$ (——) that can be achieved at different Reynolds numbers. (b) Loci of unstable poles ($\times$) and critical zeros ($\bullet$, blue; $\bullet$, red) of transfer functions $\widetilde{P}(s)$ at different Reynolds numbers. (c) Lift distributions (—— for real part and - - for imaginary part) on the cylinder at RHP zeros I and II and $Re=100$.

The control of vortex shedding using such body-mounted set-up is more challenging than control with the in-flow set-up in § 4. This is revealed by closed-loop simulations in three ways. First, controllers designed for $Re=60,80$ and $100$, although stabilising, show poorer performance than controllers designed for the in-flow set-up. Second, the optimal controller fails to stabilise the flow system if the Reynolds number is greater than 100. Third, the optimal stability margin $b_{opt}$, which is a performance indicator of the controller, decreases from 0.2537 at $Re=60$ to an extremely small value of 0.0313 at $Re=100$. This is a much more severe degradation than that seen for the in-flow set-up for which the value decreased from 0.3952 at $Re=60$ to 0.2495 at $Re=100$.

A more detailed trend of the optimal stability margin $b_{opt}$ is depicted in figure 11(a) as a function of Reynolds number. The severe degradation of control performance is clearly shown by the reduction in $b_{opt}$ from 0.4484 at $Re=50$ to 0.0313 at $Re=100$, whereas the optimal stability margin $b_{opt}$ of the in-flow set-up changes from 0.7 to 0.2495 in the same range of Reynolds numbers.

From the perspective of control theory, the performance and robustness of an optimal controller, as quantified by the optimal stability margin $b_{opt}$, can be affected by the roots (zeros and poles) of the corresponding system, especially those near the unstable mode. In general, each zero blocks a specific input signal and each RHP zero blocks a specific input signal which is unbounded (Hoagg & Bernstein Reference Hoagg and Bernstein2007). If a RHP zero occurs at exactly the same location as an unstable pole, it blocks the unstable mode exactly and the instability cannot be measured by the sensor. The system is thus unobservable and no controller is able to stabilise the unstable mode.

Figure 11(b) shows the root loci near the unstable mode for the flow systems between $Re=50$ and $Re=110$. It can be seen from the figure that as Reynolds number increases, these poles and zeros move into the RHP, which implies stronger instability and more RHP zeros. More importantly, the RHP zero I moves closer to the unstable poles at higher Reynolds numbers, which further blocks the effect of instability and reduces the observability of the unstable mode. This trend can be better observed in table 1, which summarises the stability margin and the distances between the unstable pole and RHP zeros at each Reynolds number. With increasing Reynolds number, the number of RHP zeros is increasing and their distances to the unstable pole are decreasing. The signals from the unstable mode are thus more likely to be blocked by RHP zeros. In other words, the observability of the instability is becoming worse when more RHP zeros move closer to the unstable pole. These RHP zeros are thus problematic for control design and restrict the performance and robustness of optimal controllers. Similar restrictions are also observed from the root loci of the flow system with an in-flow set-up. However, at low Reynolds numbers, the in-flow set-up does not have RHP zeros, whereas in the system with body-mounted set-up, at least one RHP zero occurs near the unstable pole.

For the body-mounted set-up, the sensor and actuator are attached to the cylinder and no time delay should be present in the transfer function from the actuator to the sensor. In other words, none of these RHP zeros is caused by the time delay. To investigate the physical mechanisms represented by these zeros, we take the cylinder flow at $Re=100$ as an example and plot the lift distributions on the cylinder surface for two RHP zeros in figure 11(c). It is interesting to note that even though zero lift should act on the cylinder at these zeros, the force acting on the cylinder surface is non-zero. Instead, the vertical force acting on the rear surface ($-90^{\circ }<\unicode[STIX]{x1D703}<90^{\circ }$) cancels the force acting on the front surface. Such force cancellations always occur near the unstable pole which blocks the signal from the unstable mode and prevents it being measured (the force on the rear surface cancels the force on the front surface). By using an integrated quantity as the measurement, the observability of the unstable mode can either be enhanced or weakened. A previous study tried surface pressure measurements for bluff body control (Dalla Longa, Morgans & Dahan Reference Dalla Longa, Morgans and Dahan2017). In that study it was found that doing (pressure on rear top) minus (pressure on rear bottom) gave better performance than summing over the whole surface (i.e. a lift-type measurement).

Table 1.

Evolution with Reynolds number of stability margin $b_{opt}$ and distance between unstable pole $p_{u}$ and RHP zero $z_{i}$.