# Computing the sensitivity of drag and lift in flow past a circular cylinder: Time-stepping versus self-consistent analysis

Philippe Meliga

► To cite this version:

Philippe Meliga. Computing the sensitivity of drag and lift in flow past a circular cylinder: Time-stepping versus self-consistent analysis. Physical Review Fluids, 2017, 2 (7), pp.073905. 10.1103/PhysRevFluids.2.073905. [hal-01698605]

HAL Id: hal-01698605

https://hal.science/hal-01698605v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Computing the sensitivity of drag and lift in flow past a
circular cylinder : time-stepping vs. self-consistent analysis

Philippe Meliga

Aix Marseille Université, CNRS, Centrale Marseille,
M2P2 UMR 7340, 13451, Marseille, France

(Dated: March 31, 2017)
Abstract

We provide in-depth scrutiny of two methods making use of adjoint-based gradients to compute the sensitivity of drag in the two-dimensional, periodic flow past a circular cylinder (Re ≲ 189). First, the time-stepping analysis used in Meliga et al. (Phys. Fluids 2014; vol.26; 104101) that relies on classical Navier-Stokes modeling and figures out the sensitivity to any generic control force from time-dependent adjoint equations marched backwards in time. Second, a novel self-consistent approach building on the model of Mantič-Lugo et al. (Phys. Rev. Lett. 2014; vol.113; 084501) to compute semi-linear approximations of the sensitivity to the mean and fluctuating components of the force. Both approaches are applied to open-loop control by a small secondary cylinder, and allow identifying the sensitive regions without knowledge of the controlled states. The theoretical predictions obtained by time-stepping analysis reproduce well the results obtained by direct numerical simulation of the two-cylinder system. So do the predictions obtained by self-consistent analysis, which corroborates the relevance of the approach as a guideline for efficient and systematic control design in the attempt to reduce drag, even though the Reynolds number is not close to the instability threshold and the oscillation amplitude is not small. This is because, unlike simpler approaches relying on linear stability analysis to predict the main features of the flow unsteadiness, the semi-linear framework encompasses rigorously the effect of the control on the mean flow, but also on the finite-amplitude fluctuation that feeds back nonlinearly onto the mean flow via the formation of Reynolds stresses. Such results are especially promising as the self-consistent approach figures out the sensitivity from time-independent equations that can be solved iteratively, which makes it generally less computationally demanding. We ultimately discuss the extent to which relevant information can be gained from a hybrid modeling computing self-consistent sensitivities from the post-processing of DNS data. Application to alternative control objectives such as increasing the lift and alleviating the fluctuating drag and lift is also discussed.

Keywords: Self-consistent modeling, sensitivity analysis, adjoint methods; drag reduction
# I. INTRODUCTION

The seminal analysis of Strykowski & Sreenivasan [1] provides experimental evidence that a small circular cylinder positioned in the near wake of a main cylinder alters vortex shedding at Reynolds numbers $Re \sim 50–100$ closely above the first instability threshold. For specific locations of the control cylinder, the authors find a stabilization of the wake accompanied by a decrease of the shedding frequency that can go towards complete suppression of unsteadiness. Since then, similar results have been obtained by direct numerical simulation (DNS) and global stability analysis of the two-cylinder system performed at about the same Reynolds numbers [2–4]. The control cylinder also acts on the resultant force, with specific positions reportedly reducing the (time-averaged) mean drag and the fluctuating lift at larger but still moderate Reynolds numbers ranging from 100 to 3,000 [5, 6]. Experimentally, the method has proven successful up to high, turbulent Reynolds numbers of order $10^4–10^5$; see Refs. [7–11] tackling circular, square and D-shaped geometries of the main cylinder. Depending on the geometry, a maximum drag reduction by 20–30% is reported with either frequency increase or decrease, but shedding itself would not be extinguished on behalf of the large Reynolds numbers. Cadot, Thiria & Beaudoin [12] have also assessed the ability of a second control cylinder (i.e., a third cylinder) in further increasing the base pressure of the D-shaped cylinder, hence additionally reduced drag.

Strykowski & Sreenivasan [1] present their results in terms of sensitivity maps showing regions close to the main cylinder where shedding is most affected by the control cylinder. In the same vein, Refs. [7, 8, 11, 13] map global quantities (e.g., Strouhal number, mean or root mean square values of drag and lift), the cost of which rapidly becomes prohibitive since systematical experimental measurements, numerical simulations or stability analyses must be performed over large parameter spaces including chiefly the position and diameter of the control cylinder (to give a taste, the experimental Strouhal number maps documented in Parezanovič & Cadot [11] and Meliga, Cadot & Serre [13] assemble shedding frequencies measured at respectively $\sim 5,000$ and $\sim 12,000$ sampled positions). Only a limited number of positions (of about a few ten) is considered in all other aforementioned studies, which yields an undersampled estimate of the real optimal. As an illustration, Dalton, Xu & Owens [5] fix the gap distance separating the centers of the two cylinders to 1.4 diameters of the main cylinder,
then vary only the angle between the center-to-center line and the free stream direc-
tion, and report a maximum drag reduction by 33% in flow past a circular cylinder at
$Re = 100$. For the exact same case, Yildirim, Rindt & Steenhoven [6] fix the stream-
wise position of the control cylinder to 0.75 diameter of the main cylinder, then vary
only its cross-wise position and report a maximum drag reduction by only 6.5%. The
discrepancy of course arises from both studies spanning different near-wake regions,
and motivates the development of more systematical approaches relying on theoretical
analysis to map quickly the best positions for placement of the control cylinder.

The experiment of Strykowski & Sreenivasan [1] has been revisited theoretically by
Hill [14] in the limit of infinitely small control cylinders. The key steps can be summa-
rized as follows : first, the growth rate and eigenfrequency of the eigenmode responsible
for the onset of vortex shedding is determined from a global stability analysis of the
steady, uncontrolled cylinder flow. The control-induced variation of the eigenvalue then
arises from the inner product between a sensitivity function making use of the adjoint
method (mathematically representing the variational derivative of the eigenvalue to a
source of momentum) and a body force mimicking the presence of the control cylinder.
By doing so, Hill identifies flow regions where the so-modeled cylinder stabilizes the
unstable eigenmode, or decrease its eigenfrequency, in good qualitative agreement with
the experimental data of Strykowski & Sreenivasan [1]. Such an approach offers an
attractive alternative to bottleneck “trial and error” procedures, as it allows spanning
quickly all possible positions of the control cylinder without ever calculating the actu-
ally controlled states. It has led a substantial body of recent work focusing on steady
and unsteady effects modeling the presence of the control cylinder, [15–21] and is now
applied to a variety of laminar [22–25] and turbulent [13, 26, 27] flows as a mean to
gain valuable insight into the most sensitive regions for open-loop control based on the
underlying physics.

This research intends to provide in-depth scrutiny of two methods making use of
similar adjoint-based gradients to compute the sensitivity of the (mean and fluctuat-
ing) drag and lift forces acting on a circular cylinder in the periodic regime. Such a
**know-how** is an asset in aerodynamic applications, where the ability to reduce the
mean drag, enhance the mean lift, and alleviate the fluctuating draft and lift is key
to improving the performances. **We go back first to** the time-stepping approach used
by Meliga, Boujo, Pujals & Gallaire [28] to analyze the sensitivity of drag in flow past a square cylinder (see also Refs. [29, 30]). The time-stepping approach relies on DNS to calculate the uncontrolled cylinder flow, and figures out the sensitivity from time-dependent adjoint equations marched backwards in time. We shall see that this is computationally very demanding because the DNS solution must be available at all adjoint time-steps, and the adjoint simulation must run long enough for a time-periodic regime to develop and for the adjoint solution to reach statistical equilibrium. An alternative approach is thus proposed, that builds on a single harmonic approximation of the uncontrolled cylinder flow governed by the semi-linear model of Mantič-Lugo, Arratia & Gallaire [31]. The latter couples a quasi-static, parameterized approximation of the instantaneous mean flow to its leading eigenmode via the formation of Reynolds stresses. This sets-up a closed, self-consistent description of the mean flow/fluctuation interaction capturing accurately the main features of the saturated oscillation, as established from comparisons with DNS data up to $Re \sim 100$. In the same range of Reynolds numbers, Meliga, Boujo & Gallaire [32] report that a self-consistent sensitivity analysis successfully predicts the effect of the control cylinder on the limit-cycle frequency and amplitude, as computed either by self-consistent modeling or DNS of the two-cylinder system. The premise is that, unlike simpler methods unable to properly describe the mean flow modification induced by the growth of unstable disturbances and the nonlinear saturation of these disturbances, self-consistent sensitivity analysis can predict with similar accuracy the effect on drag and lift, and thereby offer a credible alternative to the time-stepping approach in so far as it relies on the iterative resolution of time-independent adjoint equations, with no need for a full history of time and space-accurate solutions.

The manuscript is organized as follows. The flow configuration is described in Sec. II. Reference results obtained by time-stepping sensitivity analysis are reported in Sec. III, together with DNS results of open-loop control by a small control cylinder and a discussion on the physical mechanisms at play. Section IV gets back to the main features of the semi-linear model, and assesses its relevance for drag and lift control. Section V revisits the effect of the control cylinder by self-consistent sensitivity analysis, and provides exhaustive comparison with DNS and semi-linear results obtained by self-consistent modeling of the two-cylinder system. Finally, Sec. VI compares the costs
of both methods and examines the extent to which meaningful theoretical predictions can be obtained by a hybrid modeling computing self-consistent sensitivities from the post-processing of DNS data.

## II. FLOW CONFIGURATION

We investigate the two-dimensional (2-D), periodic, incompressible flow past a span-wise infinite circular cylinder of diameter $d^*$. A Cartesian coordinate system is used with origin at the cylinder center and drag force (resp. lift force) positive in the stream-wise +x direction (resp. in the cross-wise +y direction). Constant density $\rho^*$ and kinematic viscosity $\nu^*$ is assumed, therefore the sole parameter for this problem is the Reynolds number $Re = u_\infty^*d^*/\nu^*$, with $u_\infty^*$ the free-stream velocity. The velocity vector is $\mathbf{u} = (u, v)$ with $u$ and $v$ the stream-wise and cross-wise components. Pressure is denoted by $p$. The flow motion in space domain $\Omega$ is governed by the Navier-Stokes equations (NSE)

$$ \partial_t \mathbf{u} + \mathbf{N}(\mathbf{u}) = \mathbf{0}, \qquad (1) $$

with no-slip condition $\mathbf{u} = \mathbf{0}$ at the surface of the cylinder, denoted by $\Gamma$. In Eq. 1, $\mathbf{N}(\mathbf{u}) = \mathbf{N}(p, \mathbf{u}) = \mathbf{u} \cdot \nabla \mathbf{u} - \nabla \cdot \boldsymbol{\sigma}(p, \mathbf{u})$ is the Navier-Stokes operator, whose dependency on the pressure is omitted to ease the reading, and $\boldsymbol{\sigma}(p, \mathbf{u}) = -p\mathbf{I} + Re^{-1}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$ is the linear stress tensor. Because of incompressibility, it is understood that all velocity fields are divergence free, therefore we will not write this condition explicitly. For complex fields defined respectively in the space domain or at the surface of the cylinder, we define the inner products

$$ \langle \mathbf{v} | \mathbf{w} \rangle_{\Omega} = \int_{\Omega} \mathbf{v} \cdot \mathbf{w} \, d\Omega = \int_{\Omega} \mathbf{v}^{H} \mathbf{w} \, d\Omega, \qquad \langle \mathbf{v} | \mathbf{w} \rangle_{\Gamma} = \oint_{\Gamma} \mathbf{v} \cdot \mathbf{w} \, d\Gamma = \oint_{\Gamma} \mathbf{v}^{H} \mathbf{w} \, d\Gamma, \qquad (2) $$

where $\cdot$ is the Hermitian dot product and $H$ denotes the conjugate transpose.

Direct numerical simulations (DNS) of the NSE are performed using the finite elements/Crank-Nicholson solver presented in Meliga et al. [28]. We use a rectangular computational domain made of 108,018 triangles (378,660 degrees of freedom), with upstream and downstream boundaries at $x = -30$ and $x = 60$, respectively, and lateral boundaries at $y \pm 25$. Classical open flow conditions are imposed on the outer boundary of the space domain, namely a uniform free-stream at the inflow, symmetric
conditions at the lateral boundaries, and an advective condition at the outflow with zero pressure at the upper-right corner of the domain. The accuracy of the force coefficients in terms of mean and root mean square (rms) values has been checked to be within 2%.

For purposes of flow control, we mimic the approach of Strykowski & Sreenivasan [1] and assess how a small, secondary circular cylinder inserted at position $x_c$ in the flow modifies the drag and lift coefficients per unit length of the main cylinder (simply termed drag and lift) defined by

$$D = 2\langle \mathbf{i} | \mathbf{T}(\mathbf{u}) \rangle_{\Gamma}, \quad L = 2\langle \mathbf{j} | \mathbf{T}(\mathbf{u}) \rangle_{\Gamma}, \quad (3)$$

where $\mathbf{i}$ (resp. $\mathbf{j}$ is the unit vector in the stream-wise (resp. cross-stream) direction, $\mathbf{T}(\mathbf{u}) = \mathbf{T}(p, \mathbf{u}) = \sigma(p, \mathbf{u}) \cdot \mathbf{n}$ is the stress vector and $\mathbf{n}$ is the unit outward vector normal to the cylinder surface. The main focus here is on reducing the mean drag $\overline{D}$, but we also discuss alternative objectives, e.g., increasing the mean lift $\overline{L}$ or alleviating the rms of drag and lift defined as $D_{rms}^2 = \overline{D'^2}$ and $L_{rms}^2 = \overline{L'^2}$, where an overline denotes a mean quantity and a prime superscript denotes the zero-mean fluctuation. For several positions of interest, we perform numerical simulations of the two-cylinder system with no-slip condition at the surface of the control cylinder, accurately discretized using 300 points, as in Ref. [28]. In the most general case, we however invoke the limit of small diameter ratios and figure out the effect of the control from the inner product between a sensitivity function and a simple model of the force $\delta \mathbf{f}$ exerted by the cylinder on the flow, namely a force opposite to the drag felt by the control cylinder in a uniform flow at the local velocity

$$\delta \mathbf{f}(\mathbf{x}) = -\frac{1}{2} \eta D_{\eta}(Re_{\eta}) |\mathbf{u}| \delta(\mathbf{x} - \mathbf{x}_{c}). \quad (4)$$

In Eq. (4), $\delta$ is the 2-D Dirac delta function, $|\mathbf{u}|$ is the norm of the velocity induced by the dot product, $\eta$ is the diameter of the control cylinder, and $D_{\eta}$ is its drag coefficient taken to depend solely on the local Reynolds number $Re_{\eta}(\mathbf{x}) = Re_{\eta}|\overline{\mathbf{u}}|(\mathbf{x})$, which is because the advection time scale in the vicinity of the control cylinder is much smaller than the vortex shedding period (roughly by a factor of the diameter ratio), so the large-scale vortices shed from the main cylinder are essentially sensitive to the mean of the control cylinder flow. The relevance of this model for the intended application is discussed extensively in Ref. [28], where it is shown to reproduce with good accuracy the effect of the true control cylinder. Given the values of $\eta$ and $Re$ considered herein,
the local Reynolds number remains below ~ 20, a range where the relationship $D_\eta(Re_\eta)$ is approximated using the three-parameter power law defined in the above reference.

### III. TIME-STEPPING SENSITIVITY ANALYSIS

For a comprehensive representation of the equations taking part in the time-stepping approach, the reader is referred to Meliga et al. [28]. Only the main concepts are reviewed in the sequel, while Sec. VI comments further on the numerical cost. The intent here is not so much to discuss the detailed sensitivity of the circular cylinder, as to give points of comparison for further analysis of the self-consistent results while discussing the main physical mechanisms at play. As long as the meaning is clear from the context, we denote by cylinder flow either the *uncontrolled* cylinder flow, i.e., the solution to the NSE, or the *controlled* cylinder flow, i.e., the solution to the NSE with body force $\delta \mathbf{f}$ as an additional right-hand side (RHS). Quantities obtained by DNS are termed *nonlinear*, while their approximation obtained by sensitivity analysis are termed *linear*. Unless specified otherwise, all results pertain to a control cylinder of diameter $\eta = 0.1$.

#### A. Time-dependent adjoint equations

The linear sensitivity of the resultant force proceeds from inner products

$$ \overline{\langle \mathbf{u}^\dagger | \delta \mathbf{f} \rangle_\Omega} = \overline{\langle \mathbf{u}^\dagger | \delta \bar{\mathbf{f}} \rangle_\Omega} + \overline{\langle \mathbf{u}^{\dagger'} | \delta \bar{\mathbf{f}}' \rangle_\Omega}, \quad (5) $$

where $\mathbf{u}^\dagger$ is the solution to the initial-value problem

$$ - \partial_t \mathbf{u}^\dagger + \mathbf{L}^\dagger(\mathbf{u}) \mathbf{u}^\dagger = \mathbf{0}, \qquad \mathbf{u}^\dagger(t=0) = \mathbf{0}, \quad (6) $$

and $\mathbf{L}^\dagger(\mathbf{u}) = -\mathbf{u} \cdot \nabla \mathbf{u}^\dagger + \mathbf{u}^\dagger \cdot \nabla \mathbf{u}^T - \sigma(-p^\dagger, \mathbf{u}^\dagger)$ is the adjoint of the linearized Navier-Stokes operator, hence the adjoint LNSE moniker. Equation (6) comes with adjoint open flow boundary conditions under the form of homogeneous boundary conditions at the inflow, symmetric conditions at the lateral boundaries and an adjoint stress-free condition at the outflow [17]. The adjoint solution is driven by inhomogeneous boundary conditions at the cylinder surface

$$ \mathbf{u}^\dagger|_\Gamma = 2\theta[\xi + (1-\xi)\frac{D'(t)}{D_{\text{rms}}} \mathbf{i} + 2(1-\theta)[\xi + (1-\xi)\frac{L'(t)}{L_{\text{rms}}} \mathbf{j}], \quad (7) $$
<table><caption>Table I. Boolean values of $\theta$ and $\xi$ and relation to the mean and rms values of drag and lift.</caption><thead><tr><th></th><th>δD&#773;</th><th>δL&#773;</th><th>δD<sub>rms</sub></th><th>δL<sub>rms</sub></th></tr></thead><tbody><tr><td>(θ, ξ)</td><td>(1, 1)</td><td>(0, 1)</td><td>(1, 0)</td><td>(0, 0)</td></tr></tbody></table>

where we use booleans $\theta$ and $\xi$ to define the quantity subjected to the sensitivity analysis (mean or rms, drag or lift, as synthesized in Tab. I), each of which requires solving a dedicated adjoint problem. The general picture is that the sensitivity of the mean drag (resp. mean lift) is driven by a steady velocity aligned on the mean drag (resp. lift) vector, while that of the rms drag (resp. rms lift) is driven by a time-dependent velocity following the zero-mean fluctuation of the drag (resp. lift) vector.

From the numerical standpoint, one key feature of the adjoint equations is the reversal in space-time directionality, which shows through the minus sign ahead of the $\partial_t u^\dagger + u \cdot \nabla u^\dagger$ term accounting for the material derivative at the cylinder flow velocity. This implies to march Eq. (6) backwards in time, which in turn requires knowledge of the history of the cylinder flow solution through the entire time-span of the adjoint simulation, as further discussed in Sec. VIA. All results presented in the sequel are obtained using the approach described in [28], i.e., solving first the NSE, writing all time steps to disk, solving the adjoint LNSE over the same time interval and with the same time step, and discarding the early/late time steps (corresponding to transients of the DNS and adjoint solutions) to compute meaningful time-averages of the sensitivity integrands.

## B. Mean drag sensitivity

We start by assessing the effect of the control cylinder on the mean cylinder drag. The map of the variations $\delta D$ shown in Fig. 1(a) at $Re = 100$ is largely reminiscent of that documented for the square cylinder at the same Reynolds number (see Fig. 7 in [28]), with drag increase close to the mean separation points, and drag reduction either upstream of the main cylinder, close to the centerline, or downstream, within and on either side of the recirculation region. Although the non-fixed separation points of the circular cylinder are classically expected to bring improved control efficiency, the
<figure><img src="image_2.png" /><figcaption>Figure 1. (a)-(b) Variations of time-averaged, mean drag induced by a control cylinder of diameter $\eta = 0.1$ whose effect is modeled by Eq. (4) - $Re = 100$. (a) Main cylinder. (b) Two-cylinder system. The grey circles mark the positions $x_c = (1, 0.7)$ and $x_c = (0, 0.6)$ for which Fig. 2 provides a detailed comparison of the linear and nonlinear results. The dash-dotted line is the time-averaged, mean recirculating streamline.</figcaption></figure>

same maximum reduction by 20% prevails in the upstream region. This is much less than the 33% reported by Dalton et al. [5], a difference that we believe is attributable to these authors using quite arbitrarily a value $\overline{D} = 1.5$ of the uncontrolled mean drag. The present value $\overline{D} = 1.336$ is lower by 12% (hence the lesser efficiency) but obviously more consistent with the reference data $\sim 1.32 - 1.33$ of the literature; see, e.g., [33–35]. The drag of the two-cylinder system varies by

$$ \delta \overline{D}_{\text{tot}} = \delta \overline{D} + 2 \langle \mathbf{i} | \delta \overline{\mathbf{f}} \rangle_{\Omega}, \qquad (8) $$

and is essentially reduced in the recirculation region (where both effects add to one another in Eq. (8) since the drag of the main cylinder decreases and the control cylinder is a source of thrust); see the map in Fig. 1(b) exhibiting a maximum reduction by 8%. This is a noticeable difference with respect to the square cylinder, for which a maximum drag reduction by 10% is achieved in the upstream region of largest sensitivity (we get back to this matter in the next paragraph).

The effect of a control cylinder placed at $x_c = (1, 0.7)$ is further illustrated in Fig. 2(a). Large, open (resp. filled) circles are used to represent the nonlinear uncontrolled (resp. controlled) drag against $Re$ in the range [$Re_c^{2D}$; $Re_c^{3D}$], where $Re_c^{2D}$ is the critical Reynolds number for the onset of 2-D instability, determined by linear
stability analysis of the steady cylinder flow, and $Re_c^{3D}$ is the critical value for the onset of 3-D instability, determined by Floquet analysis [36, 37] of the 2-D periodic flow.<fnref n="1" /> Drag decreases continuously with $Re$, by 8% at $Re = 60$, up to 15% at $Re = 180$. Similar to previous results for forced wakes [40, 41], the interpretation can be made that the control cylinder increases the formation length, i.e., the size of the region bounded by the detached shears and the eddy roll-up [42]. The instantaneous flow at the rear of the body is less curved, which weakens the pressure gradient due to the centrifugal forces and raises the base pressure. This shows in Fig. 3(a) through the damping of the Reynolds stresses and the downstream shift of their spatial structure, as in Refs. [10, 11]. These trends carry over to the position $x_c = (0, 0.6)$ documented in Fig. 2(b), where drag increases by up to 25% at $Re = 120$. For this position, the control cylinder decreases the formation length, which enhances the pressure gradient and lowers the base pressure; see the related strengthening of the Reynolds stresses and the upstream shift of their spatial structure in Fig. 3(b). The same mechanism explains the drag reduction achieved inside the recirculation region and in the upstream region of largest sensitivity. In this last case, however, the induced adverse pressure gradient acts by lowering the upfront pressure, but leaves the base pressure essentially unchanged, which results in the controlled and uncontrolled Reynolds stresses having almost identical amplitude and position (not shown here). Compared to the square cylinder, the pressure gradient is softened by the bluff circular shape, which results in the velocity being larger by 10% in the in-between region. This inflates the drag of the control cylinder by 5%, and explains the loss of efficiency when it comes to the drag of the two-cylinder system.

For both positions, the linear approximations computed from Eq. (5) and superimposed in Fig. 2 as the small circles exhibit good agreement with the nonlinear data, which gives confidence that sensitivity analysis can be used as a systematic path to guide the best positions for placement of a control cylinder. There does exist a small discrepancy between the linear and nonlinear values, that we ascribe to our cylinder force model being excessively simple for a not-that-small diameter $\eta = 0.1$; see the

<footnote n="1">¹ All threshold values are specific to a given position of the control cylinder, since the control modifies also the growth rate of flow disturbances. For the same reason, we find the 3-D transition to be led by the classical A-mode [38] for the uncontrolled case and for the controlled case with control cylinder at $x_c = (0, 0.6)$, but by the so-called C-mode for the case with control cylinder at $x_c = (1, 0.7)$. This mode is too weak to compete with the A mode and B modes under natural conditions, but shows</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2. (a)-(b) Mean drag against $Re$ for control by a cylinder of diameter $\eta = 0.1$ at (a) $\mathbf{x}_c = (1, 0.7)$ and (b) $\mathbf{x}_c = (0, 0.6)$. Large open (resp. filled) symbols denote nonlinear values of the uncontrolled (resp. controlled) drag. Small symbols denote linear values obtained by time-stepping sensitivity analysis, with dark and light grey shades to stress the contributions of the mean and fluctuating components of the model force (4). (c)-(d) Same as (a)-(b) for a control cylinder of diameter $\eta = 0.02$.</figcaption></figure>

improved agreement reported in Figs. 2(c)-(d) using a smaller diameter $\eta = 0.02$. The grey shade used to represent the contribution of the individual inner products $\langle \mathbf{u}^\dagger | \delta \mathbf{f} \rangle_\Omega$ (dark shade) and $\langle \mathbf{u}^{\dagger'} | \delta \mathbf{f}^{\prime} \rangle_\Omega$ (light shade) suggests that the mean and fluctuating components of the control force contribute equally to the drag reduction achieved at $\mathbf{x}_c = (1, 0.7)$, while the drag increase at $\mathbf{x}_c = (0, 0.6)$ is triggered by the mean force only. The general picture is that the control cylinder acts mainly through the mean component of the force, except inside the recirculation region where the observed drag reduction is driven by the fluctuating component. There is a region of overlap on either
<figure><img src="image_4.png" /><figcaption>Figure 3. (a) x-component of the Reynolds stresses divergence for the uncontrolled cylinder flow (lower half plane) and for the flow controlled by a cylinder of diameter $\eta = 0.1$ at $\mathbf{x}_c = (1, 0.7)$ (upper half plane). (b) Same as (a) for a control cylinder at $\mathbf{x}_c = (0, 0.6)$ - $Re = 100$.</figcaption></figure>

side of the recirculation, where both components of the force are of the same order of magnitude. Similar disjointness in the sensitive regions has been reported for the limit-cycle frequency and amplitude of the circular cylinder [32].

### C. Mean lift and rms sensitivity

The mean lift variations mapped in Fig. 4(a) compare in every way to that of the square cylinder. Lift increases in the bottom shear layer spreading from the upstream stagnation point to the lower separation point, where the control essentially deflects the shear layer upwards, and raises the pressure over the lower cylinder surface. The similarity carries over to the variations of the rms, whose related maps in Fig. 4(b) unveil a decrease (resp. an increase) in the upstream centerline region and downstream, within and on either side of the recirculation region (resp. close to the separation points and in a large secondary region surrounding the recirculation). Those are the regions where the control cylinder decreases (resp. increases) the oscillation amplitude (see Fig. 14 in [32]), so it can be inferred that the mitigation of the fluctuations is correlated to a damping of the instability. Note, the map of rms lift shown in the lower half of Fig. 4(b) displays variations larger by two orders of magnitudes, which reflects the different amplitudes of the drag and lift oscillations. Given the uncontrolled values
<figure><img src="image_5.png" /><figcaption>Figure 4. (a) Variations of time-averaged, mean lift induced by a control cylinder of diameter $\eta = 0.1$. (b) Same as (a) for the rms drag (upper half) and lift (lower half) - Re = 100.</figcaption></figure>

of rms drag and lift ($D_{rms} = 0.0064$ vs. $L_{rms} = 0.23$) the achieved relative variations are actually comparable.

# IV. A SELF-CONSISTENT MODELING OF DRAG AND LIFT

## A. Model description

Mantič-Lugo et al. [31] have recently revisited the nonlinear saturation of the shedding instability using a semi-linear model that couples a quasi-static, parameterized approximation of the instantaneous mean flow to its leading eigenmode, considered a relevant guess for the first harmonic of the fluctuation. The eigenmode feeds back onto the mean flow via its Reynolds stresses, which sets up a closed, self-consistent description of the mean flow/fluctuation interaction. Several recent studies [43, 44] have breathed new life into this concept of a mean flow coupled to Reynolds stresses modeled from the averaged product of single eigenmode disturbances, that dates back to the early works of Malkus [45] and Stuart [46]. The main assumption is that the flow nonlinearity involves little production of higher harmonics, which is true of the cylinder flow [47]. Mantič-Lugo et al. [31] then proceed to show that saturation occurs when the mean flow is neutrally stable [40, 48, 49], at which point the eigenfrequency of the leading eigenmode predicts well the exact oscillation frequency, as early noticed in cylinder flows and related bluff-body wakes [50–52]. Similar results have been reported
by Thiria, Bouchet & Wesfreid [53] using a weakly non-parallel stability analysis of an instantaneous mean flow computed on-the-fly from DNS data. In contrast, the mean flow is not taken for a given in the self-consistent approach, but comes instead as an output of the model. Exhaustive comparison with uncontrolled DNS data corroborates the relevance of this approach to accurately predict the frequency and amplitude of vortex shedding, but also the spatial structure of the mean flow, the fluctuation and the Reynolds stresses, up to Reynolds number $Re \sim 100$ [31, 54]. Similar agreement has been reported recently for vortex shedding control by means of a small control cylinder [32].

In the periodic regime, the model quantities are governed by the self-consistent equations (SCE)

$$N(U_m) = A^2 \psi(\hat{u}_1), \quad (9)$$

$$(\lambda + i\omega)\hat{u}_1 + L(U_m)\hat{u}_1 = 0, \quad (10)$$

$$\lambda = 0, \quad (11)$$

$$\langle \hat{u}_1 | \hat{u}_1 \rangle_{\Omega} = 1, \quad (12)$$

with boundary conditions $U_m|_{\Gamma} = \hat{u}_1|_{\Gamma} = 0$. Here, $L(U_m)\hat{u}_1 = u \cdot \nabla U_m + \hat{u}_1 \cdot \nabla U_m - \sigma(\hat{p}_1, \hat{u}_1)$ is the linearized Navier-Stokes operator, $U_m$ is the self-consistent mean flow, $\hat{u}_1$ is the unit-norm eigenmode of growth rate $\lambda = 0$, eigenfrequency $\omega$ and real amplitude $A$, $\Re(\cdot)$ and $*$ indicate respectively the real part and the conjugate of a complex quantity, and we note $\psi(\hat{u}) = \hat{u} \cdot \nabla \hat{u}^* + \hat{u}^* \cdot \nabla u$. All quantities are termed semi-linear to stress that the mean flow is forced nonlinearly by the Reynolds stresses of the eigenmode, while the eigenmode is solution of a linear disturbance equation (not to be confused with the semi-linear terminology used to classify partial differential equations). For given Reynolds number, the analysis consists in determining the value of $A$ yielding neutral stability, which is achieved relaxing condition (11) and increasing the amplitude up to the point where $\lambda = 0$. Numerically, we resort on the finite elements solver presented in Meliga et al. [32] to solve iteratively the SCE with a combination of Newton and shift-and-invert Arnoldi algorithms, as further discussed in Sec. VIB. We take all solutions to be determined to a sufficient accuracy when $|\lambda| \le 10^{-6}$.
We augment here the above SCE with the additional disturbance equation

$$A^2[2(\lambda + i\omega)\hat{u}_2 + \mathbf{L}(\mathbf{U}_m)\hat{u}_2] = -A^2\zeta(\hat{u}_1), \quad (13)$$

with boundary condition $\hat{u}_2|_\Gamma = 0$, where we note $\zeta(\mathbf{u}) = \mathbf{u} \cdot \nabla \mathbf{u}$ the compact form of the nonlinear advection operator. This is the equation for the second harmonic $\hat{u}_2$, taken with no loss of generality with amplitude $A^2$. The latter is entirely slaved to its forcing by the self-interaction of the eigenmode at frequency $2\omega$, since the Reynolds stresses of $\hat{u}_2$ does not feed back onto the mean flow, nor does its nonlinear interaction with $\hat{u}_1$ onto the eigenmode. Numerically, we thus solve the SCE and the second harmonic equation sequentially. The necessity to include higher harmonics shows through the fact that the uncontrolled rms drag is otherwise trivially zero because of the reflected symmetry about the centerline. This is not true if the control cylinder is offset from the centerline, but relying on the sole eigenmode can then be insufficient to predict accurate rms values, which is a little counterintuitive given that the semi-linear model has been said to overlook all harmonic contributions. Strictly speaking, the working assumption is that the coupling with the higher harmonics has little effect on the mean flow/fluctuation interaction, not that the higher harmonics are zero. Therefore, there is no contradiction in ascribing part or all of the drag and lift fluctuations to higher harmonics, as long as they remain slaved to their forcing by the lower harmonics.

## B. Relevance for aerodynamic drag and lift control

Since the stress vector is linear in the flow quantities, the semi-linear mean drag and lift are simply the drag and lift of the self-consistent mean flow

$$D_m = 2\langle \mathbf{i} | \mathbf{T}(\mathbf{U}_m) \rangle_\Gamma, \quad L_m = 2\langle \mathbf{j} | \mathbf{T}(\mathbf{U}_m) \rangle_\Gamma. \quad (14)$$

The semi-linear rms drag and lift follow after some elementary algebra as

$$\hat{D}_{\text{rms}} = \sqrt{2} (A^2 |D_1|^2 + A^4 |D_2|^2)^{1/2}, \quad \hat{L}_{\text{rms}} = \sqrt{2} (A^2 |L_1|^2 + A^4 |L_2|^2)^{1/2}, \quad (15)$$

where we note $D_k = 2\langle \mathbf{i} | \mathbf{T}(\hat{u}_k) \rangle_\Gamma$ and $L_k = 2\langle \mathbf{j} | \mathbf{T}(\hat{u}_k) \rangle_\Gamma$ the contributions of $\hat{u}_1$ and $\hat{u}_2$ to the fluctuating drag and lift. We show in Fig. 5 the results obtained for the uncontrolled cylinder flow and for the controlled cases considered in Sec. III B, with
<figure><img src="image_6.png" /><figcaption>Figure 5. (a) Mean drag, (b) mean lift, (c) rms drag and (d) rms lift against $Re$: semi-linear (blue diamonds) vs. nonlinear results (red circles). Open symbols pertain to the uncontrolled cylinder flow. Dark (resp. light) filled symbols pertain to the flow controlled by a cylinder of diameter $\eta = 0.1$ at $\mathbf{x}_c = (1, 0.7)$ (resp. at $\mathbf{x}_c = (0, 0.6)$).</figcaption></figure>

control cylinder either at $\mathbf{x}_c = (1, 0.7)$ or at $\mathbf{x}_c = (0, 0.6)$. All semi-linear values shown as the blue diamonds fall onto their nonlinear counterparts (red circles), which proves that self-consistent modeling captures accurately the mean and fluctuating forces acting on the cylinder. The reported agreement advocates slaving the dynamics of the second harmonic (which should not come as a surprise, otherwise the model would fail to predict the main features of the saturated cylinder flow). Moreover we have checked all semi-linear values to remain identical (to the fourth decimal) retaining up to the fourth harmonic, which is consistent with the results of Protas & Wesfreid [56, 57] showing that (at least in the uncontrolled case) the spectra of the drag and lift fluctuations unveil significant contributions from the first two harmonics only. In contrast, retaining
the sole eigenmode miscalculates the rms drag by up to 25% for the controlled case at $x_c = (0, 0.6)$.

## V. SELF-CONSISTENT SENSITIVITY ANALYSIS

We now use self-consistent modeling to revisit the sensitivity analysis performed in Sec. III. All related quantities are simply termed linear. As long as the meaning is clear from the context, we denote by cylinder flow either the *uncontrolled* cylinder flow, i.e., the solution to the SCE, or the *controlled* cylinder flow, i.e., the solution to the same SCE with body force

$$ \delta \mathbf{F}_m(\mathbf{x}) = -\frac{1}{2} \eta D_\eta |\mathbf{U}_m| \mathbf{U}_m \delta(\mathbf{x}-\mathbf{x}_c), \quad \delta \hat{\mathbf{f}}_1(\mathbf{x}) = -\frac{A}{2} \eta D_\eta \left( |\mathbf{U}_m| \hat{\mathbf{u}}_1 + \frac{\mathbf{U}_m \cdot \hat{\mathbf{u}}_1}{|\mathbf{U}_m|} \mathbf{U}_m \right) \delta(\mathbf{x}-\mathbf{x}_c), \quad (16) $$

as additional RHSs in Eqs. (9)-(10), where $\delta \mathbf{F}_m$ (resp. $\delta \hat{\mathbf{f}}_1$) is the self-consistent approximation of the mean force (resp. of the fundamental component of the fluctuating force) exerted by the cylinder on the flow. Rigorously speaking, the second harmonic $\delta \hat{\mathbf{f}}_2$ of the fluctuating force should also be taken into consideration to compute the sensitivity of the rms drag and lift (not that of the mean values, that do not depend on $\hat{\mathbf{u}}_2$). The above first-order force model however retains a sufficient degree of complexity with respect to the intended purpose, in the sense that we have found the corrections to the rms sensitivity maps documented in the sequel to be negligible. For all that, this amounts not to assuming zero control-induced perturbation of the second harmonic, rather to slaving the controlled harmonic to its forcing by the controlled eigenmode.

### A. Self-consistent adjoint equations

In a manner similar to what has been presented in Sec. III, we show in Appendix A that all relevant variations proceed from inner products

$$ \langle \mathbf{U}_m^\dagger | \delta \mathbf{F}_m \rangle_\Omega + 2\Re\{\langle \hat{\mathbf{u}}_1^\dagger | \delta \hat{\mathbf{f}}_1 \rangle_\Omega\}, \quad (17) $$
where $\mathbf{U}_m^\dagger$, $\hat{\mathbf{u}}_1^\dagger$ and $\hat{\mathbf{u}}_2^\dagger$ are solutions to the adjoint self-consistent equations

$$ \mathbf{L}^\dagger(\mathbf{U}_m)\mathbf{U}_m^\dagger = -2\Re\{A\phi^\dagger(\hat{\mathbf{u}}_1, \hat{\mathbf{u}}_1^*) + A^2\phi^\dagger(\hat{\mathbf{u}}_2, \hat{\mathbf{u}}_2^*)\}, \quad (18) $$

$$ A[(\lambda - i\omega)\hat{\mathbf{u}}_1^\dagger + \mathbf{L}^\dagger(\mathbf{U}_m)\hat{\mathbf{u}}_1^\dagger] = -A^2\phi^\dagger(\hat{\mathbf{u}}_1, \mathbf{U}_m^\dagger) - A^2\phi^\dagger(\hat{\mathbf{u}}_1^*, \hat{\mathbf{u}}_2^\dagger), \quad (19) $$

$$ A^2[2(\lambda - i\omega)\hat{\mathbf{u}}_2^\dagger + \mathbf{L}^\dagger(\mathbf{U}_m)\hat{\mathbf{u}}_2^\dagger] = \mathbf{0}, \quad (20) $$

$$ \langle \mathbf{U}_m^\dagger | A^2 \psi(\hat{\mathbf{u}}_1) \rangle_\Omega = (1 - \xi) \left[ \theta \frac{A^2 |D_1|^2 + 2A^4 |D_2|^2}{\hat{D}_{\text{rms}}} + (1 - \theta) \frac{A^2 |L_1|^2 + 2A^4 |L_2|^2}{\hat{L}_{\text{rms}}} \right], \quad (21) $$

$$ \alpha^\dagger = -2A\langle\hat{\mathbf{u}}_1^\dagger | \hat{\mathbf{u}}_1\rangle_\Omega - 4A^2\langle\hat{\mathbf{u}}_2^\dagger | \hat{\mathbf{u}}_2\rangle_\Omega, \quad (22) $$

with inhomogeneous boundary conditions

$$ \mathbf{U}_m^\dagger|_\Gamma = 2\xi[\theta\mathbf{i} + (1 - \theta)\mathbf{j}], \quad (23) $$

$$ \hat{\mathbf{u}}_1^\dagger|_\Gamma = 2A(1 - \xi)[\theta \frac{D_1}{\hat{D}_{\text{rms}}} \mathbf{i} + (1 - \theta) \frac{L_1}{\hat{L}_{\text{rms}}} \mathbf{j}], \quad (24) $$

$$ \hat{\mathbf{u}}_2^\dagger|_\Gamma = 2A^2(1 - \xi)[\theta \frac{D_2}{\hat{D}_{\text{rms}}} \mathbf{i} + (1 - \theta) \frac{L_2}{\hat{L}_{\text{rms}}} \mathbf{j}], \quad (25) $$

and we note $\phi^\dagger(\mathbf{u}, \mathbf{v}) = -\mathbf{u} \cdot \nabla \mathbf{v} + \mathbf{v} \cdot \nabla \mathbf{u}^T$ the adjoint of the linearized advection operator $\phi(\mathbf{u}, \mathbf{v}) = \mathbf{u} \cdot \nabla \mathbf{v} + \mathbf{v} \cdot \nabla \mathbf{u}$. The above SCE comprise of three equations (18)-(20) for the adjoint mean flow $\mathbf{U}_m^\dagger$ and for the first two harmonics $\hat{\mathbf{u}}_1^\dagger$ and $\hat{\mathbf{u}}_2^\dagger$ of the adjoint fluctuation. The first harmonic $\hat{\mathbf{u}}_1^\dagger$ does not reduce to the leading adjoint eigenmode of the self-consistent mean flow because Eqs. (18)-(19) are coupled and must be solved simultaneously, together with compatibility condition (21) and normalization condition (22), whose role is to guarantee the marginal stability of the controlled solution through the adjoint scalar parameter $\alpha^\dagger$. Equation (20) is decoupled from the rest of the problem to reflect the slaving of the second harmonic, and is thus meaningful only to the rms problem (otherwise $\xi = 1$, Eq. (25) reduces to $\hat{\mathbf{u}}_2^\dagger|_\Gamma = \mathbf{0}$ and $\hat{\mathbf{u}}_2^\dagger = \mathbf{0}$ in the absence of a stirring mechanism). The reader is referred to Appendix A for further deepening regarding the derivation of these adjoint SCE, that differ from those derived in [32] for the limit-cycle frequency and amplitude by the adjoint second harmonic equation and by the specific compatibility condition and boundary conditions. The key point to notice from a numerical perspective is that all equations are independent of time, which allows resorting to the iterative algorithm detailed in Appendix B. We further elaborate on the related cost in Sec. VI.
<figure><img src="image_7.png" /><figcaption>Figure 6. (a)-(b) Variations of self-consistent mean drag induced by a control cylinder of diameter $\eta = 0.1$ whose effect is modeled by Eq. (16) - $Re = 100$. (a) Main cylinder. (b) Two-cylinder system. The grey circles mark the positions $x_c = (1, 0.7)$ and $x_c = (0, 0.6)$ for which Fig. 7 provides a detailed comparison of the linear and semi-linear results. The dash-dotted line is the self-consistent, mean recirculating streamline.</figcaption></figure>

## B. Mean drag sensitivity

We revisit first the effect of the control cylinder on the mean drag at $Re = 100$. The self-consistent variations $\delta D_m$ mapped in Fig. 6(a) reproduce accurately those obtained by the time-stepping method and documented in Fig. 1(a), with drag increase close to the separation points and drag reduction upstream of the cylinder, close to the centerline, or downstream, within and on either side of the recirculation. The magnitude of the variations is also well predicted, for instance we recover the same maximum reduction by 20% in the upstream region. The total variations

$$ \delta D_{m,tot} = \delta D_m + 2\langle \mathbf{i} | \delta \mathbf{F}_m \rangle_{\Omega}, \qquad (26) $$

mapped in Fig. 6(b) display the same agreement, although there exist small discrepancies with respect to the exact map shown in Fig. 1(b) close to the hyperbolic re attachment point. This is because the velocity drops rapidly to zero in this region, which results in all harmonics being of the same order of magnitude as the fundamental and strongly interacting with one another. The number of harmonics needed to achieve a proper approximation of the fluctuation is thus arbitrarily large, while the self-consistent fluctuations keeps being modeled from the fundamental only (not even
<figure><img src="image_8.png" /><figcaption>Figure 7. (a)-(b) Mean drag against *Re* for control by a cylinder of diameter *η* = 0.1 at (a) **x**<sub>c</sub> = (1, 0.7) and (b) **x**<sub>c</sub> = (0, 0.6). Large, open (resp. filled) diamonds denote semi-linear values of the uncontrolled (resp. controlled) drag. Small symbols denote linear values obtained by self-consistent sensitivity analysis, with dark and light grey shades to stress the contributions of the mean and fluctuating components of the model force (16). Nonlinear and linear values obtained by the time-stepping method are reported from Fig. 2 as the various red circles. (c)-(d) Same as (a)-(b) for control by a cylinder of diameter *η* = 0.02.</figcaption></figure>

the second harmonic, since the latter is slaved to its forcing by the leading eigenmode).

The detailed effect of control cylinders at **x**_c = (1, 0.7) and **x**_c = (0, 0.6) is reckoned in Fig. 7. The self-consistent sensitivity predictions (represented by the small blue diamonds) are seen to approximate well the semi-linear values (large blue diamonds) but also the nonlinear ones (large red circles) and related sensitivity results obtained by the time-stepping method (small red circles), which corroborates the ability of self-consistent sensitivity analysis to provide both qualitative and quantitative
predictions, and thereby to guide accurately the best positions for placement of the control cylinder. The individual contributions of the various mean and fluctuating components of the control force are especially well predicted, as evidenced by the dark and light grey shades stressing the contribution of the inner products $\langle \mathbf{U}_m^\dagger | \delta \mathbf{F}_m \rangle_\Omega$ and $2\Re\{\langle \hat{\mathbf{u}}_1^\dagger | \delta \hat{\mathbf{f}}_1 \rangle_\Omega\}$. The only noticeable difference is for a control cylinder at $\mathbf{x}_c = (0, 0.6)$ for $Re > 150$, as we find the achieved drag reduction variation to be driven essentially by the fluctuating component of the force, while the variation computed by the time-stepping approach features equal contributions from the mean and fluctuating components. This is because subtle distortions in the sensitivity regions as $Re$ increases make the control cylinder drift from a region where the mean force yields a large drag reduction, to an in-between region where it has almost no effect (hence the fading of the dark grey shade in Fig. 7(b)), but keeps yielding a large drag reduction if shifted down a little to follow the deformation of the sensitivity region (not shown here). The contribution of the mean force is thus increasingly difficult to predict because the error committed is typically of the same order of magnitude as the effect to be predicted, hence the small discrepancy with respect to the semi-linear results, even for $\eta = 0.02$. However, the self-consistent analysis also slightly overestimates the drag reduction caused by the fluctuating force. Both errors somehow compensate to make the overall variation almost identical to that predicted by the time-stepping method.

### C. Mean lift and rms sensitivity

The results obtained by self-consistent sensitivity analysis of the mean lift and of the rms drag and lift exhibit the same consistency. We retrieve from the map in Fig. 8(a) that the control cylinder increases the mean lift at $Re = 100$ if positioned in the bottom shear region. The similarity with the exact results reported in Fig. 4(a) is striking, which is because the lift sensitivity is large only in flow regions where the magnitude of the Reynolds stresses is small. The variations of rms drag and lift mapped in Fig. 8(b) are also in good agreement with those in Fig. 4(b), although the same typical miscalculations are visible close to the hyperbolic reattachment point. As has been said at the beginning of this section, we have checked the variations $2\Re\{\langle \hat{\mathbf{u}}_2^\dagger | \delta \hat{\mathbf{f}}_2 \rangle_\Omega\}$ pertaining to the second harmonic of the fluctuating force to yield little to no modifications of the above rms maps. This suggests that the control-induced
<figure><img src="image_9.png" /><figcaption>Figure 8. (a) Variations of self-consistent mean lift induced by a control cylinder of diameter $\eta = 0.1$. (b) Variations of the self-consistent rms drag (upper half) and lift (lower half) - $Re = 100$.</figcaption></figure>

modification of the second harmonic can be slaved to that of the leading eigenmode, consistently with the fact that the augmented model with slaved harmonic is relevant to predict the rms drag and lift of the controlled cases considered herein.

# VI. DISCUSSION

## A. Cost of the time-stepping sensitivity analysis

The key steps for assessing the sensitivity of the resultant force by the time-stepping approach are summarized as follows: (i) march the NSE forward in time to compute the uncontrolled cylinder flow $\mathbf{u}$, (ii) march the adjoint LNSE backwards in time to compute the adjoint solution $\mathbf{u}^{\dagger}$, (iii) average the sensitivity integrand $\langle \mathbf{u}^{\dagger} | \delta \mathbf{f} \rangle_{\Omega}$ over an appropriate time-span. The approach essentially offers the advantage it is meant to offer, i.e., it allows predicting the effect of the control with no need to compute the controlled state, which represents tremendous time saving compared to covering of entire parameter spaces by DNS. It also has the advantage of correctness, in the sense that the smallness of the control amplitude is the only underlying assumption. The method especially supports arbitrarily large amplitude of oscillations, and can thus be performed at Reynolds numbers not close to the instability threshold, unlike many sensitivity analyses dedicated to assessing the effect of the control on the flow stability
[17, 21, 58, 59].

Its major drawback lies in the fact that the adjoint LNSE use the entire time history
of the cylinder flow solutions, which turns to be very computationally demanding.
Because of the time-reversal feature, marching the adjoint equations in time indeed
requires solving first the NSE and to make the DNS solution available at every adjoint
time-step, which is inconvenient and very costly in terms of storage. To give a taste, we
use here 70Gb to write to disk 600 time units of DNS solutions, plus 60Gb to store the
adjoint solutions relevant to the mean drag problem and another 200Gb to store the
adjoint solutions relevant to the mean lift and rms problems, which is not mandatory
but recommended to achieve time-efficient post-processing. In comparison, it would
take up to 5Tb to handle only the mean drag problem in a three-dimensional (3-D) case
distributing 40 points in the span-wise direction, which is beyond feasible limits. It
is of course possible to save memory using specific check-point schemes [60] consisting
in storing the solution only at sampled intervals (the so-called check-points) and in
subsequently recomputing on-the-fly from the closest check-point, but this comes at
the price of higher CPU-time.

B. Cost of the self-consistent sensitivity analysis

The key steps for assessing the sensitivity of the resultant force by the self-consistent
approach are similarly summarized as follows: (i) solve the SCE to compute the
self-consistent cylinder flow {U_m, û_1, û_2}, (ii) solve the adjoint SCE to compute the
self-consistent adjoint solutions {U_m†, û_1†, û_2†}, (iii) compute the sensitivity integrands
⟨U_m† | δF_m⟩_Ω and ⟨û_1† | δf_1⟩_Ω. From a numerical perspective, the time-independence of
the equations is a significant advantage, as it allows resorting to iterative techniques.
The corrections made at each iteration must be under-relaxed to maintain the stability
of the iterative solution, but even so, the need to store entire time history of solutions
is eliminated, hence marginal storage requirements. In terms of CPU cost, it takes up
to a few hundreds of iterations to compute the adjoint solutions to a precision of 10⁻¹²
(in L² norm), which can be achieved within a coupled of hours on a regular sequential
workstation. This is because all operators involved in Eqs. (18)-(19) are linear in the
unknowns, therefore, it suffice to perform the necessary matrix inversions beforehand
(2 for the mean drag and lift, 3 for the rms values), following which the resolution
involves merely the repeated evaluation of matrix-vector products.

The cost of the approach is roughly that of computing the self-consistent cylinder flow, which has been found to depend highly on the Reynolds number [54]. The difficulty lies in the fact that the operators in Eqs. (9)-(10) are nonlinear in the unknowns, which requires to perform repeated matrix inversions. This is worsened by the fact that the model amplitude A itself is an unknown, hence the need to repeatedly solve the system for increasing values of A up to the point where $\lambda = 0$. In practice, the cost of solving the SCE to a precision of $10^{-12}$ (in $L^2$ norm of the various solutions) is less than that of performing the DNS on the same mesh by more than 80% at $Re = 60$ (using relaxation factors of $\sim 0.8$) and roughly 40% at $Re = 100$ (using relaxation factors of $\sim 0.3$). For values $Re > 130$, the number of iterations increases rapidly because the strong nonlinearities at play considerably slow down the convergence rate of iterative methods, as discussed for instance in [61]. At $Re = 150$, we use a relaxation factor $\sim 0.1$, which makes the cost of solving once the SCE about 3 times larger than at $Re = 100$. However, the saturation amplitude being also larger, the cost of converging the self-consistent solutions to the same precision is about 5 times larger, which ends up being twice the cost of performing the DNS. For $Re = 180$, this figure increases to 5 times, and no convergence could be achieved beyond $Re = 190$, as if the 3-D instability about to show-up somehow left a trace on the self-consistent solution. For all that, there is true added value in resorting to self-consistent sensitivity analysis : the cost remains lower than that of the time-stepping approach over a significant range of Reynolds numbers $Re < 120 - 130$, and while it is faster to compute the cylinder flow by DNS for $Re > 130$, the advantage of the self-consistent approach regarding the resolution of the adjoint equations remains unmatched.

From the accuracy standpoint, the fact that self-consistent modeling provides accurate sensitivity predictions constitutes a promising achievement compared to simpler approaches relying on classical, linear stability theory to predict the main features of the flow unsteadiness. As an example, a so-called mean flow approach is proposed in Meliga et al. [28] to compute approximated sensitivity results from knowledge of only the mean cylinder flow, the cost of which is marginal since one solves only once a steady equation for the adjoint mean flow

$$L^{\dagger}(U_{m}^{\dagger}) = 0, \qquad (27)$$
<figure><img src="image_10.png" /><figcaption>Figure 9. Variations of time-averaged mean drag induced by a control cylinder of diameter $\eta = 0.1$, approximated in the mean flow approach solving Eqs. (27)-(28) - Re = 100.</figcaption></figure>

with boundary condition

$$U_m^\dagger|_\Gamma = 2\theta \mathbf{i} + 2(1-\theta)\mathbf{j}. \qquad (28)$$

The approximation is harsh, however, because it amounts to set the harmonic of the adjoint fluctuation to zero (as can be seen comparing to Eqs. (18)-(23)) and thus inherently predicts zero rms sensitivity. For the mean drag and lift, the accuracy of the predictions is questionable, as seen from Fig. 9 showing the so-computed variations of the cylinder drag at Re = 100. This is because the approach fails to encompass the nonlinear interaction between the mean flow and its fluctuation. Actually, it does allow the control to modify both, but it does not allow the modification of the fluctuation to feed back onto the mean, which is troublesome because the frequency and amplitude of the oscillation, and thus the distribution of drag and lift, are precisely driven by the flow response to the Reynolds stresses of the fluctuation. The upstream sensitive region yielding the maximum drag reduction is well predicted in Fig. 9 (consistently with the fact that the upstream amplitude of the Reynolds stresses is negligible), but the method is clearly exposed in the wake as it misses on the drag reduction occurring in the recirculation region (where the amplitude of the Reynolds stresses peaks). Moreover, the sensitive region associated with drag reduction in the shear layers comes with indisputably different location and spatial extension; see Figs. 1(a) and 6(a) for comparison.
### C. Hybrid DNS/self-consistent sensitivity analysis

We examine now the extent to which relevant predictions can be obtained from a hybrid approach combining DNS to compute the cylinder flow solution, and self-consistent modeling to compute the sensitivities, which we believe is an interesting compromise in terms of cost. We assume that only a few snapshots of the cylinder flow (not the entire time history) are available from a DNS, together with an accurate time-averaged mean flow $\overline{U}$ used to approximate the self-consistent mean flow $U_m$. We proceed to compute the leading eigenmode and its eigenfrequency from classical mean flow stability analysis, and solve the eigenvalue problem

$$ (\lambda_{dns} + i\omega_{dns})\hat{\mathbf{u}}_{1dns} + \mathbf{L}(\overline{\mathbf{U}})\hat{\mathbf{u}}_{1dns} = \mathbf{0}, \quad (29) $$

with boundary condition $\hat{\mathbf{u}}_{1dns}|_{\Gamma} = \mathbf{0}$ and normalization condition $\langle \hat{\mathbf{u}}_{1dns} | \hat{\mathbf{u}}_{1dns} \rangle_{\Omega} = 1$. Since the cylinder mean flow is close to neutral stability [52], the obtained growth rate $\lambda_{dns}$ is automatically close to its self-consistent counterpart $\lambda = 0$. We leave aside the second harmonic because the focus is only on the mean drag sensitivity, but it is straightforward to compute the related approximation $\hat{\mathbf{u}}_{2dns}$ from knowledge of the eigenmode. For the model amplitude, we compute first the exact Reynolds stresses divergence as

$$ -\overline{\psi(\mathbf{u}')} = \mathbf{N}(\overline{\mathbf{U}}), \quad (30) $$

where the RHS can be evaluated from knowledge of the sole mean flow, then the Reynolds stresses divergence of the leading eigenmode as $\psi(\hat{\mathbf{u}}_{1dns})$, and set the amplitude to the least-square value

$$ A_{dns} = \frac{\langle \psi(\hat{\mathbf{u}}_{1dns}) | \psi(\mathbf{u}') \rangle_{\Omega}}{\langle \psi(\mathbf{u}') | \psi(\mathbf{u}') \rangle_{\Omega}}, \quad (31) $$

minimizing the difference between the exact DNS solution and its expected self-consistent approximation. Finally, we solve an adjoint system formally identical to (18)-(21), only all approximations determined above have substituted for the exact, self-consistent quantities.

We show in the upper half of Fig. 10(a) a map of the mean drag variations $\delta D_m$ computed at $Re = 100$, using for simplicity the same body forces as in (16) to model the presence of the control cylinder. The results exhibit reasonable agreement with those obtained by self-consistent and time-stepping analysis, reproduced in the lower half of
<figure><img src="image_11.png" /><figcaption>Figure 10. (a) Variations of self-consistent mean drag induced by a control cylinder of diameter $\eta = 0.1$. The variations in the upper half stem from a hybrid approach relying on adjoint equations formally identical to (18)-(21), only extrapolated DNS approximations have substituted for the exact, self-consistent solutions. The variations in the lower half are those obtained by full self-consistent sensitivity analysis, reported for comparison from Fig. 6(a). (b) Variations of time-averaged, mean drag obtained by the time-stepping method, reported for comparison from Fig. 1(a) - $Re = 100$.</figcaption></figure>

Fig. 10(a) and in Fig. 10(b), respectively. Nonetheless, there exist some discrepancies in the magnitudes of sensitivity, as the hybrid method overestimates by almost 50% the drag variations achieved in the downstream shear regions. The overall agreement is encouraging, although the miscalculations tend to worsen as the Reynolds number increases; see Fig. 11 providing similar comparisons at $Re = 150$ and $Re = 180$. Further evidence is provided in Fig. 12 showing the detailed effect of the control cylinder at $x_c = (1, 0.7)$ and $x_c = (0, 0.6)$. While the drag increase at $x_c = (0, 0.6)$ and $Re < 130$ is predicted within 10%, the drag reduction at $x_c = (1, 0.7)$ is twice as large as expected. The effect of the control is increasingly mispredicted as the Reynolds number increases. At $Re > 170$, even the sign of the variation is erroneous, but this traces back to subtle distortions of the sensitivity regions similar to those discussed in Sec. V B, for which the approximations of hybrid modeling serve as a mere magnifier (see Fig. 11 showing that the control cylinder drifts from a region where it yields a large drag reduction, to the periphery of a region where it yields a large drag increase, but keeps yielding a large drag reduction if shifted down a little to follow the deformation of the sensitivity
<figure><img src="image_12.png" /><figcaption>Figure 11. Same as Fig. 10 for (a)-(b) Re = 150 and (c)-(d) Re = 180. Approximated results obtained by the mean flow approach are shown in the close-ups.</figcaption></figure>

region).

Note, there is no unique way to determine the model amplitude. It is relevant as well to rely on the mean oscillation amplitude defined as the rms of the fluctuation energy, to give

$$A_{dns}^2 = \frac{\overline{\langle \mathbf{u}' | \mathbf{u}' \rangle_{\Omega}}}{\langle \hat{\mathbf{u}}_{1dns} | \hat{\mathbf{u}}_{1dns} \rangle_{\Omega}}. \qquad (32)$$

This is actually simpler than going through the Reynolds stresses, but we insist that the Reynolds stresses can be evaluated a posteriori from the sole mean flow, while the oscillation amplitude needs to be computed on-the-fly. Anyhow, the sensitivity results are barely impacted, as indicated by the yellow symbols in Fig. 12. This suggests that the miscalculated variations result from subtle differences between the exact, self-consistent solutions {$U_m, \hat{u}_1$} and their approximations {$\overline{U}, \hat{u}_{1dns}$}, not from an inaccuracy of the model amplitude. These differences are discussed extensively by
Mantič-Lugo, Arratia & Gallaire [54], who report that the mean flows differ essentially far downstream, in a region where the flow dynamics consists mainly in the advection and diffusion of the shed vortices. This yields a velocity deficit that is more concentrated in the time-averaged mean flow $\overline{U}$ than in its self-consistent counterpart $U_m$, hence a stronger shear that yields larger amplitudes of its eigenmode $\hat{u}_{1ds}$, but also of the adjoint solutions $\{U_m^\dagger, \hat{u}_1^\dagger\}$, and thus of the sensitivities. In a sense, this warns agains the temptation to retain only the mathematical aspect of the self-consistent equations while ignoring the underlying physics : on the one hand, the subtle failure of the self-consistent mean flow in correctly predicting the far field distribution of the time-averaged mean flow is necessary to ensure a better approximation of the vortex shedding structure based on a single eigenmode. On the other hand, the structure of the self-consistent eigenmode is the only one whose Reynolds stresses force the self-consistent mean flow in a manner such that it generates exactly the self-consistent eigenmode. Even so, we believe the main result to keep in mind is that the approach somehow predicts relevant information regarding the location of the sensitive regions. As such, it constitutes a significant achievement over the mean flow approach, whose related maps are shown in the close-ups of Figs. 11(a) and (c).

## VII. CONCLUSION

We use various adjoint-based methods to compute sensitivity maps for the (mean and fluctuating) drag and lift of a circular cylinder, a classical time-stepping analysis that figures out the exact sensitivity to generic control force from adjoint equations marched backwards in time, and a novel self-consistent analysis that builds on the model of Mantič-Lugo et al. [31] to compute semi-linear approximations of the sensitivity to the mean and fluctuating components of the force. Both approaches are applied to open-loop control by a small secondary cylinder, whose presence in the flow is modeled by a reacting force localized at the same location where the control cylinder is placed, equal and opposite to the anticipated drag.

Consistently with sensitivity results reported for flow past a square cylinder [28], the time-stepping approach provides a relevant and systematic guideline to map the best positions for placement of the control cylinder, as established from comparisons with dedicated DNS of the two-cylinder system. It is however very demanding be-
<figure><img src="image_13.png" /><figcaption>Figure 12. (a)-(b) Drag against $Re$ for control by a cylinder of diameter $\eta = 0.1$ at (a) $\mathbf{x}_c = (1, 0.7)$ and (b) $\mathbf{x}_c = (0, 0.6)$. Large, open (resp. filled) diamonds denote semi-linear values of the uncontrolled (resp. controlled) drag. Small symbols denote linear values obtained by hybrid self-consistent/DNS modeling, as obtained determining the model amplitude from the Reynolds stresses of the DNS (blue symbols) or its mean oscillation amplitude (yellow symbols). Nonlinear values are reported from Fig. 2 as the various red circles.</figcaption></figure>

cause a DNS solution of the uncontrolled cylinder flow must be available at all adjoint time-steps, and the adjoint simulation must be run long enough for a time-periodic state to show up and for the adjoint solution to reach statistical equilibrium, which yields considerable CPU and storage costs. The same agreement carries over to the predictions obtained by self-consistent analysis, which constitutes a promising achievement compared to simpler but harsher methods overlooking the nonlinear coupling between the mean and fluctuating components of the flow modification induced by the control. Moreover, the self-consistent approach figures out all sensitivities from time-independent adjoint equations, which allows resorting to efficient, iterative algorithms and reduces the cost to that of computing the self-consistent, uncontrolled cylinder flow. This is faster than performing the DNS for Reynolds numbers up to $Re = 120 - 130$. For $Re > 130$, the DNS is faster, but the flaws of the time-stepping approach regarding the resolution of the adjoint equations specifics remains, so the choice of a particular method is essentially a question of balancing constraints.

Finally, we show that valuable information regarding the sensitivity regions can be gained from a hybrid method combining DNS to compute the cylinder flow, and self-
consistent analysis to compute the sensitivity, although the magnitude of sensitivity in
the shear layers is systematically overestimated because of seemingly minor differences
between the self-consistent and time-averaged mean flows. This stresses the need
to devote further efforts to improve the numerical resolution of the self-consistent
equations (for instance using predictor-corrector techniques to anticipate the direction
of the corrections), as the related sensitivity has proven to be a valuable asset to
efficiently and accurately control unsteady flows, even though the Reynolds numbers
is not close to the instability threshold and the oscillation amplitude is not small. This
is especially true given that the relevance of such self-consistent models, coupling a
perturbation equation linearized around a mean flow forced in return by the Reynolds
stresses of the perturbation, goes well beyond periodic flows. Theoretical analysis
of turbulent flows is another important scope, for which the time-stepping approach
is bound to fail because the high sensitivity to initial conditions yields exponentially
diverging solutions if the length of the adjoint simulation exceeds the predictability time
scale [29]. We believe the time independence of the adjoint self-consistent equations
makes it a promising theory to build on in this regard, but further research is required to
assess how the small-scale turbulence feeds back on the mean flow/coherent structures
interaction, whether it be under the form of an eddy viscosity or of a stochastic forcing
in the linearized fluctuation equation. There is hope, however, as the latter approach
has recently led to savvy interpretations of how streamwise rolls and streaks self-sustain
in linearly stable wall-bounded shear flows [43, 62, 63].

ACKNOWLEDGMENTS

This work is supported by the «Investissements d'Avenir» French Government
program, managed by the French National Research Agency (ANR) through the
A*MIDEX grant (ANR-11-IDEX-0001-02) and the LABEX MEC project (ANR-11-
LABX-0092).
Appendix A: Derivation of the self-consistent sensitivity functions

We start from the augmented, forced SCE

$$
\mathbf{N}(\mathbf{U}_m) = -A^2\psi(\hat{\mathbf{u}}_1) + \mathbf{F}_m, \quad (A1)
$$

$$
A[(\lambda + i\omega)\hat{\mathbf{u}}_1 + \mathbf{L}(\mathbf{U}_m)\hat{\mathbf{u}}_1] = \hat{\mathbf{f}}_1, \quad (A2)
$$

$$
A^2[2(\lambda + i\omega)\hat{\mathbf{u}}_2 + \mathbf{L}(\mathbf{U}_m)\hat{\mathbf{u}}_2] = -A^2\zeta(\hat{\mathbf{u}}_1), \quad (A3)
$$

$$
\lambda = 0, \tag{A4}
$$

$$
\langle \hat{\mathbf{u}}_1 | \hat{\mathbf{u}}_1 \rangle_{\Omega} = 1, \qquad (A5)
$$

with homogeneous conditions $\mathbf{U}_m|_{\Gamma} = \hat{\mathbf{u}}_1|_{\Gamma} = \hat{\mathbf{u}}_2|_{\Gamma} = \mathbf{0}$ at the cylinder surface. We resort here to Lagrangian optimization, using the self-consistent solution $\mathcal{V}_d = \{\mathbf{U}_m, \hat{\mathbf{u}}_1, \hat{\mathbf{u}}_2, \lambda, \omega, A\}$ to Eqs. (A1)-(A5) as so-called direct variables, and the body forces $\mathcal{V}_c = \{\mathbf{F}_m, \hat{\mathbf{f}}_1\}$ as control variables. We recall the definitions of the mean and rms components of the self-consistent drag

$$
D_m = 2\langle \mathbf{i} | \mathbf{T}_m \rangle_{\Gamma}, \qquad \hat{D}_{\text{rms}}^2 = 2 \sum_{k=1}^{2} A^{2k} |D_k|^2, \qquad D_k = 2\langle \mathbf{i} | \mathbf{T}_k \rangle_{\Gamma}, \quad (A6)
$$

where we use $\mathbf{T}_m = \mathbf{T}(\mathbf{U}_m)$ and $\mathbf{T}_k = \mathbf{T}(\hat{\mathbf{u}}_k)$ for compact notation. In the sequel, we note similarly $\delta D_m = 2\langle \mathbf{i} | \delta \mathbf{T}_m \rangle_\Gamma = 2\langle \mathbf{i} | \mathbf{T}(\delta \mathbf{U}_m) \rangle_\Gamma$ where $\delta \mathbf{U}_m$ is the control-induced, linear perturbation to the mean flow, and $\delta D_k = 2\langle \mathbf{i} | \delta \mathbf{T}_k \rangle_\Gamma = 2\langle \mathbf{i} | \mathbf{T}(\delta \hat{\mathbf{u}}_k) \rangle_\Gamma$ where $\delta \hat{\mathbf{u}}_{1,2}$ are the related perturbations to the leading eigenmode and to the slaved harmonic. All lift-related quantities by substituting $\mathbf{j}$ (resp. $L_k$) for $\mathbf{i}$ (resp. $D_k$).

It is convenient to use the compound variable

$$
H = \xi[\theta D_m + (1 - \theta)L_m] + (1 - \xi)[\theta \hat{D}_{\text{rms}} + (1 - \theta)\hat{L}_{\text{rms}}], \quad (A7)
$$

that represents either the mean or rms value of drag and lift, depending on the value
of the Boolean doublet $(\theta, \xi)$ given in Tab. I. We seek sensitivity functions $\nabla_{\mathbf{F}_m} H$ and
$\nabla_{\hat{\mathbf{f}}_1} H$ such that

$$
\begin{align}
\delta H &= \xi[\theta \delta D_m + (1-\theta)\delta L_m] + 2(1-\xi) \sum_{k=1}^{2} A^k \Re\{[\theta \frac{D_k^*}{\hat{D}_{\text{rms}}} \delta D_k + (1-\theta) \frac{L_k^*}{\hat{D}_{\text{rms}}} \delta L_k]\}, \nonumber \\
&= \langle \nabla_{\mathbf{F}_m} H | \delta \mathbf{F} \rangle_{\Omega} + 2\Re\{\langle \nabla_{\hat{\mathbf{f}}_1} H | \delta \hat{\mathbf{f}}_1 \rangle_{\Omega}\}, \tag{A8}
\end{align}
$$
for any small modification {$\delta \mathbf{F}, \delta \hat{\mathbf{f}}_1$} of the control forces. At this stage, we introduce adjoint variables $V_a = \{\mathbf{U}_m^\dagger, \hat{\mathbf{u}}_1^\dagger, \hat{\mathbf{u}}_2^\dagger, \alpha^\dagger, \beta^\dagger\}$ and define the functional

$$ \mathcal{L}(V_c, V_s, V_a) = H - \langle \mathbf{U}_m^\dagger | \mathbf{N}(\mathbf{U}_m) + A^2 \psi(\hat{\mathbf{u}}_1) - \mathbf{F}_m \rangle_\Omega \\ -2\Re\{\langle \hat{\mathbf{u}}_1^\dagger | A[(\lambda + i\omega)\hat{\mathbf{u}}_1 + \mathbf{L}(\mathbf{U}_m)\hat{\mathbf{u}}_1] - \hat{\mathbf{f}}_1 \rangle_\Omega\} \\ -2\Re\{\langle \hat{\mathbf{u}}_2^\dagger | A^2[2(\lambda + i\omega)\hat{\mathbf{u}}_2 + \mathbf{L}(\mathbf{U}_m)\hat{\mathbf{u}}_2 + \zeta(\hat{\mathbf{u}}_1)] \rangle_\Omega\} \\ -\alpha^\dagger\lambda - \beta^\dagger(1 - \langle \hat{\mathbf{u}}_1 | \hat{\mathbf{u}}_1 \rangle_\Omega), \quad (A9) $$

that is real since the first and second harmonics $\hat{\mathbf{u}}_1^\dagger$ and $\hat{\mathbf{u}}_2^\dagger$ of the adjoint fluctuation are complex, while the adjoint mean flow $\mathbf{U}_m^\dagger$ and the adjoint scalar parameters $\alpha^\dagger$ (ensuring neutral stability) and $\beta^\dagger$ (ensuring unit norm of the eigenmode) are real. The functional reduces to $H$ if the direct variables are solutions to the SCE (A1)-(A5), in which case all differentials (in the sense of Gateau) with respect to the adjoint variables are zero. Assuming the differentials with respect to the direct variables to be similarly zero, the total variation is given by

$$ d\mathcal{L} = \langle \mathbf{U}_m^\dagger | \delta \mathbf{F}_m \rangle_\Omega + 2\Re\{\langle \hat{\mathbf{u}}_1^\dagger | \delta \hat{\mathbf{f}}_1 \rangle_\Omega\} = \delta H, \quad (A10) $$

and it follows from (A8) that $\nabla_{\mathbf{F}_m} H = \mathbf{U}_m^\dagger$ and $\nabla_{\hat{\mathbf{f}}_1} H = \hat{\mathbf{u}}_1^\dagger$. The stationarity with respect to the direct variables is achieved using the divergence theorem to integrate by parts the RHS in (A9). The adjoint equations

$$ \mathbf{L}^\dagger(\mathbf{U}_m)\mathbf{U}_m^\dagger = -2\Re\{A\phi^\dagger(\hat{\mathbf{u}}_1, \hat{\mathbf{u}}_1^{\dagger*}) + A^2\phi^\dagger(\hat{\mathbf{u}}_2, \hat{\mathbf{u}}_2^{\dagger*})\}, \quad (A11) $$

$$ A[(\lambda - i\omega)\hat{\mathbf{u}}_1^\dagger + \mathbf{L}^\dagger(\mathbf{U}_m)\hat{\mathbf{u}}_1^\dagger] = -A^2\phi^\dagger(\hat{\mathbf{u}}_1, \mathbf{U}_m^\dagger) - A^2\phi^\dagger(\hat{\mathbf{u}}_1^*, \hat{\mathbf{u}}_2^\dagger) + \beta^\dagger\hat{\mathbf{u}}_1, \quad (A12) $$

$$ A^2[2(\lambda - i\omega)\hat{\mathbf{u}}_2^\dagger + \mathbf{L}^\dagger(\mathbf{U}_m)\hat{\mathbf{u}}_2^\dagger] = 0, \quad (A13) $$

$$ \langle \mathbf{U}_m^\dagger | A^2 \psi(\hat{\mathbf{u}}_1) \rangle_\Omega + \Re\{\langle \hat{\mathbf{u}}_1^\dagger | \hat{\mathbf{f}}_1 \rangle_\Omega\} = (1 - \xi) \sum_{k=1}^{2} k A^{2k} [\theta \frac{|D_k|^2}{\hat{D}_{\text{rms}}} + (1 - \theta) \frac{|L_k|^2}{\hat{L}_{\text{rms}}}], \quad (A14) $$

$$ \alpha^\dagger = -2A\langle\hat{\mathbf{u}}_1^\dagger | \hat{\mathbf{u}}_1\rangle_\Omega - 4A^2\langle\hat{\mathbf{u}}_2^\dagger | \hat{\mathbf{u}}_2\rangle_\Omega. \quad (A15) $$

are obtained canceling all terms on $\Omega$. Canceling the bilinear concomitant on the outer boundary $\partial\Omega \setminus \Gamma$ yields homogeneous conditions at the inflow, symmetric conditions at the lateral boundaries, and an adjoint stress-free outflow condition [17]. The only terms to survive are those

$$ \langle 2\xi[\theta\mathbf{i} + (1 - \theta)\mathbf{j}] - \mathbf{U}_m^\dagger | \delta\mathbf{T}_m \rangle_\Gamma + 2(1 - \xi) \sum_{k=1}^{2} \Re\{\langle 2A^k[\theta \frac{D_k}{\hat{D}_{\text{rms}}} \mathbf{i} + (1 - \theta) \frac{L_k}{\hat{L}_{\text{rms}}} \mathbf{j}] - \hat{\mathbf{u}}_k^\dagger | A^k \delta\mathbf{T}_k \rangle_\Gamma \}, \quad (A16) $$
on the cylinder surface, whose cancellation yields inhomogeneous conditions

$$ \mathbf{U}_m^\dagger|_\Gamma = 2\xi[\theta\mathbf{i} + (1-\theta)\mathbf{j}], \qquad \hat{\mathbf{u}}_k^\dagger|_\Gamma = 2(1-\xi)A^k[\theta\frac{D_k}{\hat{D}_{\text{rms}}}\mathbf{i} + (1-\theta)\frac{L_k}{\hat{L}_{\text{rms}}}\mathbf{j}]. \quad (\text{A17}) $$

As has been said in the body article, Eq. (20) is fully decoupled from the rest of the problem to reflect the slaving of the second harmonic. It follows that $\hat{\mathbf{u}}_2^\dagger$ is trivially zero if $\xi = 1$, i.e., if the sensitivity analysis is applied to either the mean drag or lift. The value of $\beta^\dagger$ is determined taking the inner product of Eq. (A12) with $\hat{\mathbf{u}}_1$ and integrating by parts, to give

$$ \langle \mathbf{U}_m^\dagger | A^2 \psi(\hat{\mathbf{u}}_1) \rangle_\Omega + \langle \hat{\mathbf{u}}_1^\dagger | \hat{\mathbf{f}}_1 \rangle_\Omega = -2A^2 \langle \hat{\mathbf{u}}_2^\dagger | \zeta(\hat{\mathbf{u}}_1) \rangle_\Omega + \langle \hat{\mathbf{u}}_1^\dagger | A\mathbf{T}_1 \rangle_\Gamma + \beta^\dagger. \quad (\text{A18}) $$

Substituting the LHS of (A3) for $\zeta(\hat{\mathbf{u}}_1)$ and further integrating by parts, the RHS of (A18) becomes

$$ 2\langle \underbrace{A^2[2(\lambda - i\omega)\hat{\mathbf{u}}_2^\dagger + \mathbf{L}^\dagger(\mathbf{U}_m)\hat{\mathbf{u}}_2^\dagger]}_{=0} | \hat{\mathbf{u}}_2 \rangle_\Omega + \langle \hat{\mathbf{u}}_1^\dagger | A\mathbf{T}_1 \rangle_\Gamma + 2\langle \hat{\mathbf{u}}_2^\dagger | A^2\mathbf{T}_2 \rangle_\Gamma + \beta^\dagger. \quad (\text{A19}) $$

Using conditions (A17) to calculate the terms over $\Gamma$ and retaining the real part yields

$$ \langle \mathbf{U}_m^\dagger | A^2 \psi(\hat{\mathbf{u}}_1) \rangle_\Omega + \Re\{\langle \hat{\mathbf{u}}_1^\dagger | \hat{\mathbf{f}}_1 \rangle_\Omega\} = (1-\xi) \sum_{k=1}^2 k A^{2k} \left[ \theta \frac{|D_k|^2}{\hat{D}_{\text{rms}}} + (1-\theta) \frac{|L_k|^2}{\hat{L}_{\text{rms}}} \right] + \beta^\dagger, \quad (\text{A20}) $$

and ultimately $\beta^\dagger = 0$ comparing to (A14). Since we investigate the sensitivity of the uncontrolled limit cycle ($\mathbf{F}_m = \hat{\mathbf{f}}_1 = \mathbf{0}$), the adjoint SCE reduce to

$$ \mathbf{L}^\dagger(\mathbf{U}_m)\mathbf{U}_m^\dagger = -2\Re\{A\phi^\dagger(\hat{\mathbf{u}}_1, \hat{\mathbf{u}}_1^*) + A^2\phi^\dagger(\hat{\mathbf{u}}_2, \hat{\mathbf{u}}_2^*)\}, \quad (\text{A21}) $$

$$ A[(\lambda - i\omega)\hat{\mathbf{u}}_1^\dagger + \mathbf{L}^\dagger(\mathbf{U}_m)\hat{\mathbf{u}}_1^\dagger] = -A^2\phi^\dagger(\hat{\mathbf{u}}_1, \mathbf{U}_m^\dagger) - A^2\phi^\dagger(\hat{\mathbf{u}}_1^*, \hat{\mathbf{u}}_2^\dagger), \quad (\text{A22}) $$

$$ A^2[2(\lambda - i\omega)\hat{\mathbf{u}}_2^\dagger + \mathbf{L}^\dagger(\mathbf{U}_m)\hat{\mathbf{u}}_2^\dagger] = \mathbf{0}, \quad (\text{A23}) $$

$$ \langle \mathbf{U}_m^\dagger | A^2 \psi(\hat{\mathbf{u}}_1) \rangle_\Omega = (1 - \xi) \sum_{k=1}^2 k A^{2k} \left[ \theta \frac{|D_k|^2}{\hat{D}_{\text{rms}}} + (1 - \theta) \frac{|L_k|^2}{\hat{L}_{\text{rms}}} \right], \quad (\text{A24}) $$

$$ \alpha^\dagger = -2A\langle\hat{\mathbf{u}}_1^\dagger | \hat{\mathbf{u}}_1\rangle_\Omega - 4A^2\langle\hat{\mathbf{u}}_2^\dagger | \hat{\mathbf{u}}_2\rangle_\Omega. \quad (\text{A25}) $$

## Appendix B: Resolution of the adjoint self-consistent equations

Prior to solving the above adjoint SCE, we use the Arnoldi method to compute the leading adjoint eigenmode of the self-consistent mean flow, such that

$$ (\lambda - i\omega)\hat{\mathbf{u}}_{10}^\dagger + \mathbf{L}^\dagger(\mathbf{U}_m)\hat{\mathbf{u}}_{10}^\dagger = \mathbf{0}, \quad (\text{B1}) $$
with boundary conditions $\hat{\mathbf{u}}_{10}^{\dagger}|_{\Gamma} = \mathbf{0}$, normalized to either $\langle \hat{\mathbf{u}}_{10r}^{\dagger} | \hat{\mathbf{u}}_{1} \rangle_{\Omega} = 1/2A$ or $\langle \hat{\mathbf{u}}_{10i} | \hat{\mathbf{u}}_{1} \rangle_{\Omega} = i/2A$. We also compute the adjoint mean solutions to

$$L^{\dagger}(U_m)U_{0r,i}^{\dagger} = -2 \Re\{A\phi^{\dagger}(\hat{\mathbf{u}}_1, \hat{\mathbf{u}}_{10r,i}^{*})\}, \quad (B2)$$

with homogeneous conditions $U_{0r,i}^{\dagger}|_{\Gamma} = \mathbf{0}$, as well as the inner products $\chi_{r,i} = \langle U_{0r,i}^{\dagger} | A^2\psi(\hat{\mathbf{u}}_1) \rangle_{\Omega}$ with the Reynolds stresses of the eigenmode. While $\hat{\mathbf{u}}_{10i}^{\dagger} = -i\hat{\mathbf{u}}_{10r}^{\dagger}$, there is no such relation between $U_{0r}^{\dagger}$ and $U_{0i}^{\dagger}$ because only the real part of the interaction with $\hat{\mathbf{u}}_1$ shows up in the RHS of Eq. (B2). **If the analysis relates to the rms, we also compute the adjoint second harmonic $\hat{\mathbf{u}}_2^{\dagger}$ and the inner product $\gamma_2 = \gamma_{2r} + i\gamma_{2i} = \langle \hat{\mathbf{u}}_2^{\dagger} | \hat{\mathbf{u}}_2 \rangle_{\Omega}/4A^2$ (otherwise $\hat{\mathbf{u}}_2^{\dagger}$, and thus $\gamma_2$, are trivially zero).**

Since $\hat{\mathbf{u}}_1^{\dagger}$ is solution to the forced eigenvalue problem (A22), it decomposes into the sum of homogeneous and particular solutions according to

$$\hat{\mathbf{u}}_1^{\dagger} = -(\alpha^{\dagger} + \gamma_{2r})\hat{\mathbf{u}}_{10r}^{\dagger} - \gamma_{2i}\hat{\mathbf{u}}_{10i}^{\dagger} + \hat{\mathbf{u}}_{1p}^{\dagger}, \quad (B3)$$

where the particular solution $\hat{\mathbf{u}}_{1p}^{\dagger}$ is by construction such that $\langle \hat{\mathbf{u}}_{1p}^{\dagger} | \hat{\mathbf{u}}_1 \rangle_{\Omega} = 0$ and the coefficients ahead of the homogeneous adjoint solutions deduce from (A25). The adjoint mean flow equation is linear in both $U_m^{\dagger}$ and $\hat{\mathbf{u}}_1^{\dagger}$, therefore its general solution is given in general form by

$$U_m^{\dagger} = -(\alpha^{\dagger} + \gamma_{2r})U_{0r}^{\dagger} - \gamma_{2i}U_{0i}^{\dagger} + U_p^{\dagger}, \quad (B4)$$

where $U_p^{\dagger}$ is solution to

$$L^{\dagger}(U_m)U_p^{\dagger} = -2 \Re\{A\phi^{\dagger}(\hat{\mathbf{u}}_1, \mathbf{u}_{1p}^{*}) + A^2\phi^{\dagger}(\hat{\mathbf{u}}_2, \hat{\mathbf{u}}_2^{*})\}, \quad (B5)$$

with boundary condition

$$U_p^{\dagger}|_{\Gamma} = 2\xi[\theta\mathbf{i} + (1 - \theta)\mathbf{j}]. \quad (B6)$$

Starting from a certain guess for $U_m^{\dagger}$ (for instance, $U_{0r}^{\dagger}$), we solve first Eq. (A22) and deduce the particular component from the obtained numerical solution as

$$\hat{\mathbf{u}}_{1p}^{\dagger} = \hat{\mathbf{u}}_{1}^{\dagger} - 2A\langle\hat{\mathbf{u}}_{1}^{\dagger} | \hat{\mathbf{u}}_{1} \rangle_{\Omega}\hat{\mathbf{u}}_{10r}^{\dagger}. \quad (B7)$$

This allows solving Eq. (B5) for $U_p^{\dagger}$ and computing the inner product $\chi_p = \langle U_p^{\dagger} | A^2\psi(\hat{\mathbf{u}}_1) \rangle_{\Omega}$, at which point the value of $\alpha^{\dagger}$ proceeds from the compatibility conditions (A24) to give

$$\chi_r \alpha^{\dagger} = \chi_p - \gamma_{2r} \chi_r + \gamma_{2i} \chi_i - (1 - \xi) \sum_{k=1}^{2} k A^{2k} \left[ \theta \frac{|D_k|^2}{\hat{D}_{\text{rms}}} + (1 - \theta) \frac{|L_k|^2}{\hat{L}_{\text{rms}}} \right]. \quad (B8)$$
A new guess for $\mathbf{U}_m^\dagger$ is built from Eq. (B4), and the process repeats until the difference between two consecutive iterations is less than $10^{-12}$ in $L^2$ norm (for both $\mathbf{U}_m^\dagger$ and $\hat{\mathbf{u}}_1^\dagger$).

[1] P.J. Strykowski and K.R. Sreenivasan, “On the formation and suppression of vortex ‘shedding’ at low Reynolds numbers,” J. Fluid Mech. **218**, 71–107 (1990).

[2] H.B. Kim and K.S. Chang, “Numerical study on vortex shedding from a circular cylinder influenced by a nearby control wire,” Comput. Fluid Dyn. J. **4**, 151–164. (1995).

[3] S. Mittal and A. Raghuvanshi, “Control of vortex shedding behind circular cylinder for flows at low Reynolds numbers,” Int. J. Numer. Methods Fluids **35**, 421–447 (2001).

[4] M. Morzynski, K. Afanasiev, and F. Thiele, “Solution of the eigenvalue problems resulting from global non-parallel flow stability analysis,” Comput. Meth. Appl. Mech. Engng. **169**, 161–176 (1999).

[5] C. Dalton, Y. Xu, and J.C. Owen, “The suppression of lift on a circular cylinder due to vortex shedding at moderate Reynolds numbers,” J. Fluids Struct. **15**, 617–628 (2001).

[6] I. Yildirim, C.C.M. Rindt, and A.A. Steenhoven, “Vortex dynamics in a wire-disturbed cylinder wake,” Phys. Fluids **22**, 094101 (2010).

[7] H. Sakamoto, K. Tan, and H. Haniu, “An optimum suppression of fluid forces by controlling a shear layer separated from a square prism,” J. Fluids Eng. **113**, 183–189 (1991).

[8] H. Sakamoto and H. Haniu, “Optimum suppression of fluid forces acting on a circular cylinder,” J. Fluids Eng. **116**, 221–227 (1994).

[9] T. Igarashi, “Drag reduction of a square prism by flow control using a small rod,” J. Wind Eng. Ind. Aerodyn. **69–71**, 141–153 (1997).

[10] V. Parezanović and O. Cadot, “The impact of a local perturbation on global properties of a turbulent wake,” Phys. Fluids **21**, 071701 (2009).

[11] V. Parezanović and O. Cadot, “Experimental sensitivity analysis of the global properties of a two-dimensional turbulent wake,” J. Fluid Mech. **693**, 115–149 (2012).

[12] O. Cadot, B. Thiria, and J.-F. Beaudoin, “Passive drag control of a turbulent wake by local disturbances,” in *IUTAM Symposium on Unsteady Separated Flows and their*
*Control*, IUTAM Bookseries, Vol. 14, edited by M. Braza and K. Hourigan (Springer, 2009) pp. 529–537.

[13] P. Meliga, O. Cadot, and E. Serre, “Experimental and theoretical sensitivity analysis of turbulent flow past a square cylinder,” Flow Turbul. Combust. **97**, 987–1015 (2016).

[14] D.C. Hill, “A theoretical approach for analyzing the restabilization of wakes,” NASA technical memorandum No. 103858 (1992).

[15] F. Giannetti and P. Luchini, “Structural sensitivity of the first instability of the cylinder wake,” J. Fluid Mech. **581**, 167–197 (2007).

[16] P. Luchini, F. Giannetti, and J. Pralits, “Structural sensitivity of linear and nonlinear global modes,” AIAA 2008-4227 (2008).

[17] O. Marquet, D. Sipp, and L. Jacquin, “Sensitivity analysis and passive control of cylinder flow,” J. Fluid Mech. **615**, 221–252 (2008).

[18] O. Marquet, D. Sipp, L. Jacquin, and J.-M. Chomaz, “Multiple scale and sensitivity analysis for the passive control of the cylinder flow,” AIAA 2008-4228 (2008).

[19] P. Luchini, F. Giannetti, and J. Pralits, “Structural sensitivity of the finite-amplitude vortex shedding behind a circular cylinder.” in *IUTAM Symposium on Unsteady Separated Flows and their Control*, IUTAM Bookseries, Vol. 14, edited by M. Braza and K. Hourigan (Springer, 2009) pp. 151–160.

[20] J.O. Pralits, L. Brandt, and F. Giannetti, “Instability and sensitivity of the flow around a rotating circular cylinder,” J. Fluid Mech. **650**, 513–536 (2010).

[21] P. Meliga and J.-M. Chomaz, “Global modes in a confined impinging jet: application to heat transfer and control,” Theor. Comput. Fluid Dyn. **25**, 179–193 (2011).

[22] P. Meliga, D. Sipp, and J.-M. Chomaz, “Open-loop control of compressible afterbody flows using adjoint methods,” Phys. Fluids **22**, 054109 (2010).

[23] D. Sipp, O. Marquet, P. Meliga, and A. Barbagallo, “Dynamics and control of global instabilities in open-flows: a linearized approach,” App. Mech. Rev. **63**, 030801 (2010).

[24] O. Tammisola, “Oscillatory sensitivity patterns for global modes in wakes,” J. Fluid Mech. **701**, 251–277 (2012).

[25] P. Luchini and A. Bottaro, “Adjoint equations in stability analysis.” Annu. Rev. Fluid Mech. **46**, 493–517 (2014).

[26] P. Meliga, G. Pujals, and E. Serre, “Sensitivity of 2-D turbulent flow past a D-shaped
cylinder using global stability," Phys. Fluids **24**, 061701 (2012).

[27] C. Mettot, D. Sipp, and H. Bézard, "Quasi-laminar stability and sensitivity analyses for turbulent flows: Prediction of low-frequency unsteadiness and passive control," Phys. Fluids **26**, 045112 (2014).

[28] P. Meliga, E. Boujo, G. Pujals, and F. Gallaire, "Sensitivity of aerodynamic forces in laminar and turbulent flow past a square cylinder," Phys. Fluids **26**, 104101 (2014).

[29] Q. Wang and J.-H. Gao, "The drag-adjoint field of a circular cylinder wake at reynolds numbers 20, 100 and 500," J. Fluid Mech. **730**, 145-161 (2013).

[30] X. Mao, "Sensitivity of forces to wall transpiration in flow past an aerofoil," Proc. R. Soc. A **471**, 20150618 (2015).

[31] V. Mantič-Lugo, C. Arratia, and F. Gallaire, "Self-consistent mean flow description of the nonlinear saturation of the vortex shedding in the cylinder wake," Phys. Rev. Lett. **113**, 084501 (2014).

[32] P. Meliga, E. Boujo, and F. Gallaire, "A self-consistent formulation for the sensitivity analysis of finite amplitude vortex shedding in the cylinder wake," J. Fluid Mech. **800**, 327-357 (2016).

[33] R.D. Henderson, "Details of the drag curve near the onset of vortex shedding," Phys. Fluids **7**, 2102-2104 (1995).

[34] A. G. Kravchenko and P. Moin, "B-spline method and zonal grids for simulations of complex turbulent flows," J. Comp. Phys. **151**, 757-789 (1999).

[35] S. Mittal, "Excitation of shear layer instability in flow past a cylinder at low reynolds number," Int. J. Numer. Methods Fluids **49**, 1147-1167 (2005).

[36] B.R. Noack and H. Eckelmann, "A global stability analysis of the steady and periodic cylinder wake," J. Fluid Mech. **270**, 297-330 (1994).

[37] R.D. Henderson and D. Barkley, "Secondary instability in the wake of a circular cylinder," Phys. Fluids **8**, 1683-1685 (1996).

[38] C.H.K. Williamson, "Vortex dynamics in the cylinder wake," Annu. Rev. Fluid Mech. **28**, 477-539 (1996).

[39] H.-Q. Zhang, U. Fey, B.R. Noack, M. König, and H. Eckelmann, "On the transition of the cylinder wake," Phys. Fluids **7**, 779-794 (1995).

[40] B. Thiria and J.E. Wesfreid, "Stability properties of forced wakes," J. Fluid Mech. **579**,
137–161 (2007).

[41] B. Thiria, O. Cadot, and J.-F. Beaudoin, “Passive drag control of a blunt trailing edge cylinder,” J. Fluids Struct. **25**, 766–776 (2009).

[42] J. H. Gerrard, “The mechanics of the formation region of vortices behind bluff bodies,” J. Fluid Mech. **25**, 401–413 (1966).

[43] B.F. Farrell and P.J. Ioannou, “Dynamics of streamwise rolls and streaks in turbulent wall-bounded shear flow,” J. Fluid Mech. **708**, 149–196 (2012).

[44] J.O. Pralits, A. Bottaro, and S. Cherubini, “Weakly nonlinear optimal perturbations,” J. Fluid Mech. **785**, 135–151 (2015).

[45] W. V. R. Malkus, “Outline of a theory of turbulent shear flow,” J. Fluid Mech. **1**, 521–539 (1956).

[46] J.T. Stuart, “Non-linear stability theory,” Annu. Rev. Fluid Mech. **3**, 347–370 (1971).

[47] J. Dušek, P. Le Gal, and P. Fraunié, “A numerical and theoretical study of the first Hopf bifurcation in a cylinder wake,” J. Fluid Mech. **264**, 59–80 (1994).

[48] A. Maurel, V. Pagneux, and J.E. Wesfreid, “Mean-flow correction as non-linear saturation mechanism,” Europhys. Lett. **32**, 217–222 (1995).

[49] B.J.A. Zielinska, S. Goujon-Durand, J. Dušek, and J.E. Wesfreid, “Strongly nonlinear effect in unstable wakes,” Phys. Rev. Lett. **79**, 3893–3896 (1997).

[50] D. A. Hammond and L. G. Redekopp, “Global dynamics of symmetric and asymmetric wakes,” J. Fluid Mech. **331**, 231–260 (1997).

[51] B. Pier, “On the frequency selection of finite-amplitude vortex shedding in the cylinder wake,” J. Fluid Mech. **458**, 407–417 (2002).

[52] D. Barkley, “Linear analysis of the cylinder wake mean flow,” Europhys. Lett. **75**, 750–756 (2006).

[53] B. Thiria, G. Bouchet, and J.E. Wesfreid, “On the relation between linear stability analysis and mean flow properties in wakes,” e-print arXiv:1506.05948 [physics.flu-dyn].

[54] V. Mantič-Lugo, C. Arratia, and F. Gallaire, “A self-consistent model for the saturation dynamics of the vortex shedding around the mean flow in the unstable cylinder wake,” Phys. Fluids **27**, 074103 (2015).

[55] P. Meliga, “Harmonics generation and the mechanics of saturation in flow over an open cavity : a second-order, self-consistent description,” J. Fluid Mech. (submitted) (2016).
[56] B. Protas and J.E. Wesfreid, “Drag force in the open-loop control of the cylinder wake in the laminar regime,” Phys. Fluids **14**, 810–826 (2002).

[57] B. Protas and J.E. Wesfreid, “On the relation between the global modes and the spectra of drag and lift in periodic wake flows,” C.R. Mec. **331**, 49–54 (2003).

[58] L. Brandt, D. Sipp, J.O. Pralits, and O. Marquet, “Effect of base-flow variation in noise amplifiers: the flat-plate boundary layer,” J. Fluid Mech. **687**, 503–528 (2011).

[59] E. Boujo, U. Ehrenstein, and F. Gallaire, “Open-loop control of noise amplification in a separated boundary layer flow,” Phys. Fluids **25**, 124106 (2013).

[60] A. Griewank and A. Walther, “An implementation of checkpointing for the reverse or adjoint mode of computational differentiation,” ACM T. Math. Software **26**, 19–45 (2000).

[61] B. Fornberg, “A numerical study of steady viscous flow past a circular cylinder,” J. Fluid Mech. **98**, 819–855 (1980).

[62] S. Nazarenko, N. K. Kevlahan, and B. Dubrulle, “Nonlinear rdt theory of near-wall turbulence,” Phys. D **139**, 158–176 (2000).

[63] F. Bouchet, C. Nardini, and T. Tangarife, “Kinetic theory of jet dynamics in the stochastic barotropic and 2d navier-stokes equations,” J. Stat. Phys. **153**, 572–625 (2013).