1. Introduction
Drag reduction in flow passing a circular cylinder is a canonical problem in fluid mechanics as it is related to extensive engineering applications and incorporates rich fluid physics, including the separation of boundary layers, vortex shedding, hydrodynamic stabilities, etc. Numerous control techniques have been developed for drag reduction in flow around cylinders and, depending on the energy input, they can be broadly classified as active and passive controls (Choi, Jeon & Kim Reference Choi, Jeon and Kim2008).
In active drag reduction, spanwise-uniform forcing imposed on the cylinder surface, denoted as two-dimensional (2-D) control in this work, has been extensively investigated by testing various actuations such as blowing/suction, surface rotation, etc. For example, the blowing control has been shown to reduce drag 14 % in flow passing a circular cylinder at a blowing flow rate 0.145 (normalized by the free-stream flow rate) and a maximum relative control velocity (with respect to the free-stream velocity as will be used in the following) approximately 0.37 at Reynolds number 
$Re=90$. A temporally sinusoidal rotation of a circular cylinder at 
$Re=15\,000$ has been reported to reduce drag up to 80 % in an experimental work (Tokumaru & Dimotakis Reference Tokumaru and Dimotakis1991) and 70 % in a similar 2-D simulation work (Shiels & Leonard Reference Shiels and Leonard2001). Later Fujisawa, Tanahashi & Srinivas (Reference Fujisawa, Tanahashi and Srinivas2005) conducted experimental measurement of forces acting on a rotating circular cylinder at 
$Re=2000$ and reported drag magnification at a high rotating frequency (five times the vortex shedding frequency) and 30 % drag reduction at a low rotating frequency (vortex shedding frequency). Instead of using prescribed rotary velocity, Lu et al. (Reference Lu, Qin, Teng and Li2011) used the lift coefficient as the feedback signal to calculate the rotating velocity of a circular cylinder and achieved drag reduction over 10 % at 
$Re=200$. Most recently, Palkin et al. (Reference Palkin, Hadziabdic, Mullyadzhanov and Hanjalic2018) achieved a drag reduction up to 88 % in flow passing a circular cylinder at 
$Re=1.4 \times 10^5$ by rotating the cylinder at the optimal combination of frequency and magnitude. Rather than imposing a uniform rotation, Shukla & Arakeri (Reference Shukla and Arakeri2013) studied the optimal distribution of surface rotation control around the whole surface of a circular cylinder to minimize the power loss in the limit of 
$Re \ll 1$. The obtained optimal control was observed to reduce drag up to 77 % at a maximum relative control velocity 1.5 and 
$Re=100$.
Spanwise-dependent forcing, referred to as three-dimensional (3-D) control in the following, has been shown to be more effective than 2-D control (Park et al. Reference Park, Lee, Jeon, Hahn, Kim, Kim, Choi and Choi2006; Guercio, Cossu & Pujals Reference Guercio, Cossu and Pujals2014a). Poncet et al. (Reference Poncet, Hildebrand, Cottet and Koumoutsakos2008) found that a 2-D optimal control is no longer optimal in 3-D domains, and a spanwise-varying streamwise control is more effective owing to the generation of streamwise vortex braids in a circular cylinder flow at 
$300 \le Re \le 1000$. Kim & Choi (Reference Kim and Choi2005) applied steady blowing/suction varying sinusoidally in the spanwise direction on the upper and lower surfaces of a cylinder and observed that a symmetric (‘in-phase’) control forcing reduces drag significantly and also suppresses vortex shedding more efficiently than base bleeding. They achieved approximately 20 % drag reduction at a maximum relative control velocity 0.1 at 
$Re=100$. The mechanism of this control is attributed to the phase mismatch in the spanwise direction, similar to effects induced by spanwise geometry variation.
In passive drag reduction, the configuration of the body is modified and the control action does not require energy inputs. Such shape modifications target at the separating boundary layer or the unsteady vortex shedding in the wake for drag reduction. For example, dimples on the cylinder surface induces pairs of vortices to energize the boundary layer and subsequently reduce drag from 
$Re=4\times 10^4$ to 
$Re=3\times 10^5$ (Bearman & Harvey Reference Bearman and Harvey1993). Spanwise waviness introduced to either the front or rear stagnation face has been shown to suppress vortex shedding and reduce drag in flow around both circular and rectangular cylinders (Bearman & Owen Reference Bearman and Owen1998; Darekar & Sherwin Reference Darekar and Sherwin2001; Owen, Bearman & Szewczyk Reference Owen, Bearman and Szewczyk2001). The control mechanism revealed in these works is similar to the 3-D active control reviewed above. Splitter plates, either rigid or flexible, have been also tested owing to its simple geometry and effectiveness in drag reduction and the suppression of lift fluctuations (Bearman Reference Bearman1965; Bagheri, Mazzino & Bottaro Reference Bagheri, Mazzino and Bottaro2012). Other surface modifications, including helical strakes, bumps, small tabs and so forth, were also reported to be effective in controlling the flow around a cylinder by attenuating the vortex shedding, as summarized in a review paper (Choi et al. Reference Choi, Jeon and Kim2008).
In all these efforts, the control forcing is either prescribed or computed from an optimization process or a sensitivity analysis without iterative optimizations, which can be regarded as a suboptimal computation (Min & Choi Reference Min and Choi1999). Sensitivity analyses in fluid dynamics have been focused on the structural sensitivity of stability to feedback forcing (Giannetti & Luchini Reference Giannetti and Luchini2007; Pralits, Brandt & Giannetti Reference Pralits, Brandt and Giannetti2010), and the sensitivity of modal and non-modal instabilities with respect to base flow modifications (Brandt et al. Reference Brandt, Sipp, Pralits and Marquet2011; Lashgari et al. Reference Lashgari, Tammisola, Citro, Juniper and Brandt2014; Mao Reference Mao2015a). The calculation involves the integration of the adjoint of the linearized Navier–Stokes (NS) equation, which has been widely used in non-normality, sensitivity and receptivity studies (Chomaz Reference Chomaz2005; Schmid Reference Schmid2007; Mao, Blackburn & Sherwin Reference Mao, Blackburn and Sherwin2012; Boujo & Gallaire Reference Boujo and Gallaire2014; Luchini & Bottaro Reference Luchini and Bottaro2014; Meliga et al. Reference Meliga, Boujo, Pujals and Gallaire2014), as well as in the shape optimization and optimal control investigations (Mohammadi & Pironneau Reference Mohammadi and Pironneau2004; Jameson & Ou Reference Jameson and Ou2010; Bueno-Orovio et al. Reference Bueno-Orovio, Castro, Palacios and Zuazua2012). In sensitivity analyses of flow passing a bluff body, Hwang, Kim & Choi (Reference Hwang, Kim and Choi2013) computed the sensitivity of absolute instability with respect to spanwise wavy base flow modifications, Guercio, Cossu & Pujals (Reference Guercio, Cossu and Pujals2014b) studied the optimal spanwise periodic blow/suction that suppresses parallel wake instability, and Tammisola et al. calculated the second-order sensitivity of instability to blowing/suction and surface roughness (Tammisola et al. Reference Tammisola, Giannetti, Citro and Juniper2014; Tammisola Reference Tammisola2017). In these sensitivity investigations, the base flow is 2-D and control forcing with a prescribed spanwise wavenumber was adopted to modify the base flow. It is still unclear if these sensitivity analyses can be extended to high-Reynolds-number 3-D base (or uncontrolled) flow for effective suppression of wake instabilities and drag reduction. Also it is unknown if a superposition of spanwise waves can be more efficient in controlling the wake flow than a single wave.
This work aims at answering these questions. The active control will be modelled as active boundary perturbations (velocity perturbations on the cylinder surface) while passive control will be considered as passive boundary perturbations (surface roughness). An algorithm for the sensitivity of drag to perturbations will be established. The calculation can be regarded as the first step of an optimization and only requires one call of the NS equation and the adjoint equation. This algorithm can be regarded as a combination of those developed for sensitivity of drag to external forcing (Meliga et al. Reference Meliga, Boujo, Pujals and Gallaire2014) and sensitivity of instability to boundary perturbations (Tammisola et al. Reference Tammisola, Giannetti, Citro and Juniper2014; Tammisola Reference Tammisola2017). It will be shown that when the boundary perturbation is of the passive type, the formula of the sensitivity of force is distinctively different from that of instabilities. This algorithm is then applied to flow past a circular cylinder at 
$190\le Re \le 1000$. As the spanwise-dependent controls have been demonstrated to be more effective than spanwise-independent ones, efforts will be devoted to addressing the spanwise variation of the control forcing.
2. Algorithms
2.1. Governing equations and definitions
Assuming the flow to be incompressible and Newtonian and considering the control as the boundary perturbation on the cylinder surface, the perturbed flow is governed by the NS equations:
\begin{equation} \partial_t \hat{\boldsymbol{u}}+ \hat{\boldsymbol{u}} \boldsymbol{\cdot}\boldsymbol{\nabla} \hat{\boldsymbol{u}} +\boldsymbol{\nabla} \hat{p} -Re^{{-}1}\nabla^2 \hat{\boldsymbol{u}}=0, \quad \text{with} \ \boldsymbol{\nabla} \boldsymbol{\cdot} \hat{\boldsymbol{u}}=0, \end{equation}
where 
$\hat {\boldsymbol {u}}$ and 
$\hat {p}$ are the perturbed velocity and pressure, respectively, and 
$Re$ is the Reynolds number based on the free-stream velocity and the cylinder diameter.
The perturbed flow can be decomposed as the sum of an unperturbed flow (synonymous to base flow in sensitivity, stability and non-normality studies), which is itself a solution to the NS equations, and a perturbation-induced flow, i.e. 
$(\hat {\boldsymbol {u}},\hat {p})=(\boldsymbol {U},P)+(\boldsymbol {u},p)$, where 
$\boldsymbol {U}$, 
$\boldsymbol {u}$, 
$P$ and 
$p$ are the unperturbed velocity, perturbation velocity, unperturbed pressure and perturbation pressure, respectively. Then we substitute this decomposition into the NS equations (2.1) and linearize the convection term (assuming the perturbation 
$\boldsymbol {u}$ is small enough that the interaction of this perturbation with itself is negligible) to reach the linearized NS equation
\begin{equation} \partial_t \boldsymbol{u}+ \boldsymbol{U} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U} +\boldsymbol{\nabla} p -Re^{{-}1}\nabla^2 \boldsymbol{u}=0, \quad \text{with} \ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{equation}
where the initial condition of 
$\boldsymbol {u}$ is set to zero so as to exclude effects of initial perturbations on the development of boundary perturbations. Similar to the velocity and pressure terms, the force acting on a solid body can also be decomposed into two components, namely an unperturbed force and a perturbation-induced force, as 
$\hat {\boldsymbol {f}}=\boldsymbol {F}+\boldsymbol {f}$.
To simplify notations in the following derivations, we define scalar products
\begin{align} {(}\boldsymbol{a},\boldsymbol{b})&=\int_\varOmega \boldsymbol{a}\boldsymbol{\cdot} \boldsymbol{b} \, \mathrm{d}\varOmega , \hspace{2pc} \langle \boldsymbol{a},\boldsymbol{b} \rangle=\tau^{{-}1} \int_0^\tau\int_\varOmega \boldsymbol{a}\boldsymbol{\cdot} \boldsymbol{b} \, \mathrm{d}\varOmega \,\mathrm{d}t, \end{align}
\begin{align} {[}\boldsymbol{d},\boldsymbol{e}]&=L_z^{{-}1}\int_{\boldsymbol{B}} \boldsymbol{d}\boldsymbol{\cdot} \boldsymbol{e} \, \mathrm{d}\boldsymbol{B}, \quad \{\boldsymbol{d},\boldsymbol{e}\}=L_z^{{-}1}\tau^{{-}1}\int_0^\tau \int_{\boldsymbol{B}} \boldsymbol{d}\boldsymbol{\cdot} \boldsymbol{e} \, \mathrm{d}\boldsymbol{B}\, \mathrm{d}t, \end{align}
where 
$\varOmega$ denotes the spatial domain, 
$\boldsymbol {B}$ represents the surface of the cylinder where the boundary perturbation is imposed, 
$\tau$ is a final time, 
$L_z$ is the spanwise length of the domain, 
$\boldsymbol {a},\boldsymbol {b} \in \varOmega \times [0,\tau ]$ and 
$\boldsymbol {d},\boldsymbol {e} \in \boldsymbol {B} \times [0,\tau ]$.
We also define a boundary norm:
\begin{equation} \| \boldsymbol{X}\|=[\boldsymbol{X},\boldsymbol{X}]^{1/2}, \end{equation}
which is the square root of the square integration of 
$\boldsymbol {X}$ defined on the perturbed boundary. The magnitude of a boundary perturbation can be measured by its boundary norm.
2.2. Sensitivity of forces to active boundary perturbations
In this section, the algorithm for sensitivity of forces to active control (or active boundary perturbations) is briefly introduced. The details of derivation can be found in Mao (Reference Mao2015b).
The active boundary perturbation can be modelled as Dirichlet-type velocity boundary conditions of (2.2), denoted as 
$\boldsymbol {u}(\boldsymbol {B},t)$. The spatial and temporal dependence of this perturbation is decomposed as
\begin{equation} \boldsymbol{u}(\boldsymbol{B},t)= g(t) \boldsymbol{u}_a(\boldsymbol{B}, \omega) \cos (\omega t), \end{equation}
in order to reduce the dimension of its discretized form. Here 
$\boldsymbol {u}_a(\boldsymbol {B}, \omega )$ denotes the spatial dependence of the perturbation, 
$\omega$ represents the frequency, and 
$g(t)$ is a numerical factor defined as
\begin{equation} g(t)=(1-\mathrm{e}^{-\sigma_1 t^2}). \end{equation}
This numerical factor sets 
$\boldsymbol {u}(\boldsymbol {B},0)=0$ and therefore eliminates the incompatibility of the boundary perturbation with the zero initial condition of 
$\boldsymbol {u}$. Here 
$\sigma _1$ is a relaxation factor and 
$\sigma _1=100$ is adopted throughout this work. Clearly when the final time 
$\tau \rightarrow \infty$, 
$\boldsymbol {u}_a(\boldsymbol {B})$ represents the magnitude of the boundary perturbation. This decomposition of the spatial and temporal dependences is particularly useful in open-loop flow control problems, where the control forcing is typically steady or periodic in time.
Considering that the force acting on a solid body is commonly time-dependent, when studying the perturbation-induced force, it is sensible to adopt its (quasi) mean value, whose component in direction 
$\boldsymbol {K}$ can be written as
\begin{equation} \bar{f}_{\boldsymbol{K}}= \tau^{{-}1}\int_0^\tau \boldsymbol{f} \boldsymbol{\cdot} \boldsymbol{K} R(t)\, \mathrm{d}t =\{ p\boldsymbol{n} - Re^{{-}1}\partial_{\boldsymbol{n}}\boldsymbol{u} , \boldsymbol{K} R(t) \}, \end{equation}
where 
$\boldsymbol {f}= L_z^{-1} \int _{\boldsymbol {B}} (p\boldsymbol {n} - Re^{-1} \partial _{\boldsymbol {n}}\boldsymbol {u}) \, \mathrm {d}\boldsymbol {B}$ is the perturbation-induced force per unit spanwise length. Thus 
$\bar {f}_{\boldsymbol {K}}$ can be interpreted as the perturbation-induced time-averaged force per unit spanwise length in the 
$\boldsymbol {K}$ direction. Here 
$R(t)=1-\mathrm {e}^{-\sigma _2 (t-\tau )^2}$ is a numerical factor and 
$\sigma _2=100$ is adopted throughout this work. This factor generates a zero adjoint velocity boundary condition, which is compatible with the zero initial condition of the adjoint velocity (see appendix A). Since all the variables involved in the governing equations have been non-dimensionalized, the force coefficient is two times the force.
Considering the definition of the Gâteaux differential, the gradient of the force with respect to the active perturbation, denoted as 
$\boldsymbol {\nabla }_{\boldsymbol {u}_a}\bar {f}_{\boldsymbol {K}}$, satisfies 
$\bar {f}_{\boldsymbol {K}}= [\boldsymbol {\nabla }_{\boldsymbol {u}_a}\bar {f}_{\boldsymbol {K}}, \boldsymbol {u}_a ]$. Combining (2.6), (2.8) and (A3) from adjoint analyses in appendix A, this gradient can be expressed as
\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}} =\tau^{{-}1} \int_0^{\tau} (p^*\boldsymbol{n} -Re^{{-}1} \partial_{\boldsymbol{n}}\boldsymbol{u}^* ) g \cos (\omega t)\,\mathrm{d} t, \end{equation}
where 
$p^*$ and 
$\boldsymbol {u}^*$ are solutions of the adjoint equation (A2) with appropriate initial and boundary conditions. Such an adjoint method has been extensively used in sensitivity, receptivity and non-normality calculations and shape optimization (Barkley, Blackburn & Sherwin Reference Barkley, Blackburn and Sherwin2008; Shinohara et al. Reference Shinohara, Okuda, Ito, Nakajima and Ida2008; Jameson & Ou Reference Jameson and Ou2010). If the boundary perturbation is restricted, e.g. to the wall-normal direction, the corresponding gradient should be the product of the right-hand side of (2.9) and the unit wall normal 
$\boldsymbol {n}$. It is noted that the formula of gradient in (2.9) can be also obtained by defining a Lagrangian functional and applying the variational principle.
This gradient represents the sensitivity of force with respect to active boundary perturbations. The norm of this sensitivity quantifies the controllability of the force with respect to the boundary perturbation, and the distribution of the sensitivity is parallel with the (linearly) most effective perturbation to modify the force. Therefore, the sensitivity can be scaled as a control, whose induced force is the product of the norms of the sensitivity and the control in the linear limit (Mao Reference Mao2015b).
2.3. Sensitivity of forces with respect to passive boundary perturbations
In the present linear framework, we only consider the surface-normal roughness or surface undulation, denoted as 
$r$ in figure 1(a). The more complex roughness with both surface-normal and surface-tangential components leads to cavities on the surfaces (see figure 1b). This form of roughness cannot be modelled as surface velocity perturbations under the linear assumption and will not be considered in the following. It would be possible to study this roughness by computing the sensitivity of forces with respect to the surface geometry directly.
Figure 1. Schematic plot of (a) surface-normal roughness and (b) surface-normal and surface-tangential roughness. Here 
$r$ represents the roughness. The solid lines denote the smooth surfaces without roughness 
$\boldsymbol {B}$, dashed lines represent surfaces with roughness 
$\boldsymbol {B}_r$ and dotted lines illustrate the direction of roughness.
The magnitude of the roughness can be evaluated by the boundary norm 
$\|r\|$ (see (2.5)). As shown in figure 1, after adding roughness, the surface coordinate changes from 
$\boldsymbol {B}$ to 
$\boldsymbol {B}_r=\boldsymbol {B}+r\boldsymbol { n}$, where 
$\boldsymbol {n}$ is the unit wall-normal vector towards the body.
As discussed above, when there is no roughness, the flow is unperturbed and the boundary condition on the perturbation boundary is 
$\boldsymbol {U}(\boldsymbol {B},t)=0$. When roughness is considered, the flow is the sum of the unperturbed flow and the perturbation flow, and the boundary condition is imposed on the rough surface as 
$\hat {\boldsymbol {u}}(\boldsymbol {B}_r,t)=0$. This condition can be expanded as
\begin{equation} \hat{\boldsymbol{u}}(\boldsymbol{B}_r,t)=\hat{\boldsymbol{u}}(\boldsymbol{B},t)+r \partial_{\boldsymbol{n}}\hat{\boldsymbol{u}}(\boldsymbol{B},t) +O(r^2)=0. \end{equation}
Substitute 
$\hat {\boldsymbol {u}}= \boldsymbol {U}+\boldsymbol {u}$ into this equation and linearize with respect to 
$r$ (and the perturbation velocity which is induced by the roughness) to reach
\begin{equation} \boldsymbol{u}(\boldsymbol{B},t)={-} r \boldsymbol{\partial_{\boldsymbol{n}}}\boldsymbol{U}. \end{equation}
Therefore, the flow around the rough surface can be obtained by studying an unperturbed flow and a perturbation flow, both of which have boundary conditions defined on the smooth surface 
$\boldsymbol {B}$.
When passive perturbations are considered, the surface becomes 
$\boldsymbol {B}_r$ and the force acting on this rough surface is
\begin{equation} \hat{\boldsymbol{f}}(\boldsymbol{B}_r)= \int_{\boldsymbol{B}_r} ( \hat{p} \boldsymbol{n} - Re^{{-}1} \partial_{\boldsymbol{n}} \hat{\boldsymbol{u}} ) \, \mathrm{d}\boldsymbol{B}_r, \end{equation}
where, in the integration on the right, all the variables are defined on the rough surface. It will be shown below that all the terms in the above equation (
$\hat {p}$, 
$\boldsymbol {n}$, 
$\hat {\boldsymbol {u}}$ and 
$\mathrm {d}\boldsymbol {B}_r$) can be represented by variables defined on the smooth (unperturbed) surface, provided that the magnitude of the perturbation is small enough.
As stated in § 2.1, the pressure and velocity variables at 
$\boldsymbol {B}_r$ can be decomposed as the sum of a base variable and a perturbation variable. If the boundary perturbation or roughness is small, the base and perturbation variables at 
$\boldsymbol {B}_r$ can be linearized to reach
\begin{equation} \left. \begin{gathered} \hat{p}(\boldsymbol{B}_r)=P(\boldsymbol{B}_r) + p(\boldsymbol{B}_r) \approx P(\boldsymbol{B}) + p(\boldsymbol{B})+ r \partial_{\boldsymbol{n}} P(\boldsymbol{B}),\\ \hat{\boldsymbol{u}}(\boldsymbol{B}_r) =\boldsymbol{U}(\boldsymbol{B}_r)+ \boldsymbol{u}(\boldsymbol{B}_r) \approx \boldsymbol{U}(\boldsymbol{B}) + \boldsymbol{u}(\boldsymbol{B}) + r \partial_{\boldsymbol{n}}\boldsymbol{U}(\boldsymbol{B}). \end{gathered} \right\} \end{equation}
On the far right, the first term denotes the base flow variables on the smooth surface, the second one refers to the perturbations on the smooth surface and the third term stems from the difference of the base flow on 
$\boldsymbol {B}$ and 
$\boldsymbol {B}_r$.
The unit normal defined on the rough surface can be also expressed in terms of variables defined on the smooth surface through linearization:
\begin{equation} \boldsymbol{n}(\boldsymbol{B}_r)= \boldsymbol{n}(\boldsymbol{B}) - \boldsymbol{m} \partial_{\boldsymbol{m}} r, \end{equation}
where 
$\boldsymbol {m}$ is the tangential unit vector along the smooth surface and therefore is normal to 
$\boldsymbol {n}$. The first term on the right is the ‘base’ vector and the second one is related to the roughness.
Through similar linearization, the differentiation of the rough surface can be written as a function of the differentiation of the smooth surface:
\begin{equation} \mathrm{d}\boldsymbol{B}_r= \mathrm{d}\boldsymbol{B}+ r \boldsymbol{m} \boldsymbol{\cdot} \partial_{\boldsymbol{m}} \boldsymbol{n}\,\mathrm{d}\boldsymbol{B}. \end{equation}The first term on the right is the ‘base’ differentiation of the smooth surface, while the second term is related with the roughness.
Substituting (2.13), (2.14) and (2.15) into (2.12),
\begin{equation} \hat{\boldsymbol{f}}(\boldsymbol{B}_r)=\boldsymbol{F}(\boldsymbol{B}) + \boldsymbol{f}(\boldsymbol{B}) + \boldsymbol{f}_1(\boldsymbol{B}) + \boldsymbol{f}_2(\boldsymbol{B}) + \boldsymbol{f}_3(\boldsymbol{B}), \end{equation}where
\begin{gather} \boldsymbol{F}=\int_{\boldsymbol{B}} (P\boldsymbol{n} - Re^{{-}1} \partial_{\boldsymbol{n}}\boldsymbol{U} ) \, \mathrm{d}\boldsymbol{B}, \end{gather}
\begin{gather}\boldsymbol{f}(\boldsymbol{B})= \int_{\boldsymbol{B}} ( p\boldsymbol{n} - Re^{{-}1} \partial_{\boldsymbol{n}}\boldsymbol{u} ) \, \mathrm{d}\boldsymbol{B}, \end{gather}
\begin{gather}\boldsymbol{f}_1(\boldsymbol{B})= \int_{\boldsymbol{B}} ( r \partial_{\boldsymbol{n}} P\boldsymbol{n} - Re^{{-}1} r\partial_{\boldsymbol{n}^2}\boldsymbol{U} ) \, \mathrm{d}\boldsymbol{B}, \end{gather}
\begin{gather}\boldsymbol{f}_2(\boldsymbol{B})= \int_{\boldsymbol{B}} -P \boldsymbol{m} \partial_{\boldsymbol{m}} r \, \mathrm{d}\boldsymbol{B}, \end{gather}
\begin{gather}\boldsymbol{f}_3(\boldsymbol{B})= \int_{\boldsymbol{B}} (P\boldsymbol{n} - Re^{{-}1} \partial_{\boldsymbol{n}}\boldsymbol{U} ) (r \boldsymbol{m} \boldsymbol{\cdot} \partial_{\boldsymbol{m}} \boldsymbol{n}) \, \mathrm{d}\boldsymbol{B}. \end{gather}Comparing the forces on smooth and rough surfaces, the force induced by roughness is
\begin{equation} \hat{\boldsymbol{f}}(\boldsymbol{B}_r)-\boldsymbol{F}(\boldsymbol{B}) = \boldsymbol{f}(\boldsymbol{B}) + \boldsymbol{f}_1(\boldsymbol{B}) + \boldsymbol{f}_2(\boldsymbol{B}) + \boldsymbol{f}_3(\boldsymbol{B}). \end{equation} Then the (quasi) mean force induced by roughness in direction 
$\boldsymbol {K}$ becomes
\begin{equation} \bar{f}_{\boldsymbol{K}}= \tau^{{-}1}\int_0^\tau ( \hat{\boldsymbol{f}}(\boldsymbol{B}_r)-\boldsymbol{F}(\boldsymbol{B})) \boldsymbol{\cdot} \boldsymbol{K} R(t) \,\mathrm{d}t=\{ p\boldsymbol{n} - Re^{{-}1} \partial_{\boldsymbol{n}}\boldsymbol{u} +r \boldsymbol{H}, \boldsymbol{K} R(t) \}, \end{equation}
where 
$\boldsymbol {H}=\boldsymbol {\nabla } P- Re^{-1} \partial _{\boldsymbol {n}^2}\boldsymbol {U} - Re^{-1} \partial _{\boldsymbol {n}}\boldsymbol {U} \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {n}$. Some identities are used to derive this equation, e.g. 
$\partial _{\boldsymbol {m}} \boldsymbol {m} + \boldsymbol {n} (\boldsymbol {m} \boldsymbol {\cdot } \partial _{\boldsymbol {m}}\boldsymbol {n})=0$ and 
$\boldsymbol {m} \boldsymbol {\cdot } \partial _{\boldsymbol {m}} \boldsymbol {n} = \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {n}$.
Combining (2.11), (A3) and (2.23), the mean force induced by roughness in direction 
$\boldsymbol {K}$ can be reformulated as
\begin{equation} \bar{f}_{\boldsymbol{K}}= [ \boldsymbol{\nabla}_{r}\bar{f}_{\boldsymbol{K}},r ], \end{equation}where
\begin{equation} \boldsymbol{\nabla}_{r}\bar{f}_{\boldsymbol{K}} =\tau^{{-}1} \int_0^{\tau}[ - \partial_{\boldsymbol{n}}\boldsymbol{U} \boldsymbol{\cdot} (p^*\boldsymbol{n} -Re^{{-}1} \partial_{\boldsymbol{n}}\boldsymbol{u}^* ) + R(t) \boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{H}]\, \mathrm{d} t \end{equation}
is the gradient of the mean force with respect to the roughness 
$r$. To calculate this gradient, the unperturbed flow should be computed first through integrating the NS equations and then the adjoint variables are solved by integrating the adjoint equation (A2).
This gradient represents the sensitivity of force with respect to surface roughness. Similarly as discussed in § 2.2, it can be proved that the roughness-induced force reaches maximum 
$\|\boldsymbol {\nabla }_{r}\bar {f}_{\boldsymbol {K}}\| \boldsymbol {\cdot } \|r\|$ (or minimum 
$-\|\boldsymbol {\nabla }_{r}\bar {f}_{\boldsymbol {K}}\| \boldsymbol {\cdot } \|r\|$) when 
$r$ is in the same (or opposite) direction as 
$\boldsymbol {\nabla }_{r}\bar {f}_{\boldsymbol {K}}$. The norm of the sensitivity represents the controllability of the force with respect to surface roughness, and the shape of the sensitivity around the perturbed boundary represents the optimal surface roughness, which is most effective to modify the force within the linear limit.
3. Problem description, discretization and convergence
The cylinder is centred at 
$(x,y,z)=(0,0,0)$ in the Cartesian system, and the inflow, outflow and far-field boundaries are located at 
$x=-35$, 
$x=155$ and 
$y={\pm }50$, respectively. Each 
$x$–
$y$ plane is decomposed into 2310 spectral elements, as shown in figure 2. In each element, piecewise continuous nodal-based polynomial expansions with order 
$\mathcal {P}$ are applied. For the perturbed flow, a uniform free-stream velocity condition is imposed on the inflow and far-field boundaries, and a zero-Neumann and zero-Dirichlet condition is imposed on the outflow boundary for velocity and pressure, respectively. The non-zero-Dirichlet velocity perturbation on the surface of the cylinder can be decomposed into wall-normal, azimuthal and spanwise components, as shown in figure 2(b). This velocity perturbation is set to zero in base flow calculations. Fourier decomposition is invoked in the spanwise direction and 32 Fourier modes are calculated (Blackburn & Sherwin Reference Blackburn and Sherwin2004). Time integration is carried out using a velocity-correction scheme. These numerics are adopted to integrate both the NS equations and the adjoint equation using a well-validated numerical code which has been applied in direct numerical simulations (DNS) and hydrodynamic stability studies of vortex flow and flow around solid bodies (Mao, Blackburn & Sherwin Reference Mao, Blackburn and Sherwin2013,Reference Mao, Blackburn and Sherwin2015).
Figure 2. Mesh and boundary conditions in (a) the entire domain and (b) the subdomain around the cylinder. The three components of the control are marked in (b).
In the convergence study of the sensitivity calculation, the Reynolds number is fixed at 
$Re=190$, which is slightly above the critical condition of mode A instabilities (Williamson Reference Williamson1996). Unless otherwise stated, the spanwise length is set to 
$L_z=4$, which is long enough to accommodate mode A (Barkley & Henderson Reference Barkley and Henderson1996), as shown in figure 3. We also note that at the present Reynolds number and spanwise length, the breakdown of mode A to vortex dislocations was not observed (Williamson Reference Williamson1996).
Figure 3. Iso-surfaces of spanwise vorticity 0.4 and 
$-0.4$, coloured by streamwise velocity, in unperturbed flow around a cylinder at 
$Re=190$. The thick solid lines represent the vortex lines associated with the lower vortex sheet.
To calculate the sensitivity of drag, i.e. the force in the streamwise direction, the direction 
$\boldsymbol {K}$ is set to 
$[1,0, 0]^\mathrm {T}$. Over the range of parameters studied, the sensitivities at non-zero frequencies do not yield an effective control at finite control magnitude owing to the lock-in effect (Mao Reference Mao2015b) therefore 
$\omega =0$ is adopted throughout this work.
To evaluate the convergence of the sensitivity with respect to the numerical discretization, define
\begin{equation} c_{\mathcal{P}} ({\mathcal{P}}_1, {\mathcal{P}}_2)=\frac{[\boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}}({\mathcal{P}}_1), \boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}}({\mathcal{P}}_2)]} {\| \boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}}({\mathcal{P}}_1)\| \boldsymbol{\cdot} \| \boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}}({\mathcal{P}}_2)\|}, \end{equation}
where the sensitivity is considered as a function of the polynomial order. From definition, the deviation of 
$c_{\mathcal {P}} ({\mathcal {P}}_1, {\mathcal {P}}_2)$ from 
$1$ denotes the difference of the two sensitivities obtained at polynomial orders 
${\mathcal {P}}_1$ and 
${\mathcal {P}}_2$.
Therefore, the convergence of 
$c_{\mathcal {P}} ({\mathcal {P}}, 7)$ to 
$1$ denotes the convergence of the sensitivity with respect to discretization, where the result at 
${\mathcal {P}}=7$ is taken as the reference. It is seen from table 1(a) that at polynomial order 
${\mathcal {P}}=6$, the sensitivity has converged to seven significant figures. Similarly when doubling the resolution in the spanwise direction, the change of the sensitivity is within 
$0.01\,\%$. The convergence of DNS at 
$Re=1000$, the highest 
$Re$ considered in this work, is further tested and it is found that at 
${\mathcal {P}}=6$, the mean drag converges to two significant figures with respect to that at 
${\mathcal {P}}=7$. Therefore, 
${\mathcal {P}}=6$ will be used in all the following simulations.
Table 1. Convergence of the sensitivity of drag to active perturbations in flow past a cylinder at 
$Re=190$ and 
$\omega =0$. (a) The final time is fixed at 
$\tau =40$ and the polynomial order in the spectral element method is varied. Here 
$c_{\mathcal {P}}$ is a correlation factor defined in (3.1). (b) The polynomial order is fixed at 
$\mathcal {P}=6$ and the final time is varied. Here 
$c_{\tau }(\tau _1,\tau _2)$ as defined in (3.2) denotes the correlation of sensitivities obtained at two final times.
The convergence of the calculation with respect to the final time is also studied, as shown in table 1(b). Similarly with (3.1), a correlation factor is defined to evaluate the correlation of sensitivities obtained at two final times, 
$\tau _1$ and 
$\tau _2$:
\begin{equation} c_{\tau}(\tau_1,\tau_2) = \frac{[\boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}}(\tau_1), \boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}}(\tau_2]} {\| \boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}}(\tau_1)\| \boldsymbol{\cdot} \| \boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}}(\tau_2)\|}, \end{equation}
where the sensitivity is considered as a function of the final time. Taking the sensitivity obtained at 
$\tau =200$ as the reference, the convergence of 
$c_{\tau }(\tau ,200)$ to 
$1$ denotes the convergence of the sensitivity with respect to 
$\tau$. At 
$\tau =100$, the sensitivity converges to three significant figures (doubling 
$\tau$ results in a 0.1 % change) and this final time is used in the following investigations.
4. Active drag reduction
4.1. Three-dimensional sensitivities to active perturbations
The normalized sensitivity of drag to active perturbation is decomposed into surface-normal, azimuthal and spanwise components, as shown in figures 4(a), 4(b) and 4(c), respectively, where 
$\theta$ is the azimuthal coordinate and 
$\theta=0$ at the backward stagnation point. These three components can be interpreted as the sensitivity to wall-normal, azimuthal and spanwise active perturbations, respectively. It is noticed that the sensitivity to wall-normal perturbation is much stronger than the other two, suggesting that for this flow the drag is more sensitive to wall-normal blowing/suction than surface displacement in either the azimuthal or spanwise direction. Each of the three sensitivity components consists of multiple Fourier modes in the spanwise direction. For example, the wall-normal component is a superposition of spanwise uniform blowing and periodic blowing/suction. In the azimuthal direction, all three sensitivity components are concentrated upstream of the separation lines and peak at approximately 
$\theta ={\rm \pi} /2$, similar to the control to suppress vortex shedding (Guercio et al. Reference Guercio, Cossu and Pujals2014a; Tammisola Reference Tammisola2017). These sensitivities are not ideally symmetric (figure 4a), or anti-symmetric (figure 4b,c), due to the limited number of vortex shedding periods contained in the time interval considered. It can be expected that these numerical imperfections can be eliminated when the final time 
$\tau$ is sufficiently large.
Figure 4. Contours of the sensitivity to active perturbations normalized by the boundary norm at 
$Re=190$ and 
$\omega =0$. Panels (a,b,c) are the wall-normal, azimuthal and spanwise components of the sensitivity, respectively. The dash-dotted lines represent the mean separation lines.
The spanwise length 
$L_z$ has been reported to be a critical factor in the control of cylinder flow using surface perturbations (Guercio et al. Reference Guercio, Cossu and Pujals2014a). The sensitivity of drag to wall-normal active perturbations obtained at span length 
$L_z= 4$, 
$8$ and 
$12$ is illustrated in figure 5. Other lengths including 
$L_z=2$, 
$3$, 
$5$, 
$6$ were also tested but were not shown here as the uncontrolled flow does not present 3-D instabilities. Note figure 5(a) is the same as figure 4(a) but is compressed in the spanwise direction for comparison with figures 5(b) and 5(c). In each 
$L_z$ test, the unperturbed flow is obtained from DNS over a sufficiently long time. At 
$L_z=8$ and 
$12$, the sensitivity, as well as the uncontrolled flow (not shown for brevity), is approximately a spanwise repeat of that at 
$L_z=4$, confirming that 
$L_z=4$ is sufficient for the present sensitivity analyses. This span length will be used in the following drag reduction investigations.
Figure 5. Contours of the wall-normal component of the sensitivity normalized by the boundary norm at 
$Re=190$, 
$\omega =0$ and spanwise length 4, 8 and 12 for (a,b,c), respectively. The dash-dotted lines represent the mean separation lines.
The sensitivity to active perturbations can be scaled as an active control:
\begin{equation} \boldsymbol{u}_a= S_a \frac{\boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}}}{ \|\boldsymbol{\nabla}_{\boldsymbol{u}_a}\bar{f}_{\boldsymbol{K}}\|}, \end{equation}
where 
$S_a$ is a scale factor and its absolute value denotes the boundary norm of the control. According to the definition of the sensitivity, drag reduction can be achieved by assigning a small and negative value to 
$S_a$. Since the wall-normal and tangential components of the control are generated by different actuators, only the wall-normal component or blowing/suction is considered. The control is then obtained by scaling the sensitivity of drag to wall-normal perturbations shown in figure 4(a). The spanwise distribution of the control at various azimuthal locations is plotted in figure 6(a). A scale factor 
$S_a=-7 \times 10^{-3}$ is used, and the selection of this value will be discussed in § 4.2. The control is concentrated around the centre of the domain in the spanwise direction and 
$\theta ={\rm \pi} /2$ in the azimuthal direction. Such a localized control can be generated by a limited number of actuators and will be easier to implement than the spanwise uniform or periodic controls.
Figure 6. The active control scaled from sensitivity at 
$S_a=-7 \times 10^{-3}$, 
$\omega =0$ and 
$Re=190$. (a) Spanwise distribution at various azimuthal locations and (b) dominant spanwise modes, with 
$\beta$ denoting the spanwise wavenumber.
Through Fourier decomposition in the spanwise direction, it is found that this control is dominated by three modes with spanwise wavenumber 
$\beta =0, 1$ and 
$2$ within the domain (corresponding to non-dimensionalized wavelength 
$\infty$, 
$4$ and 
$2$), whose boundary norms are 
$2.2 \times 10^{-3}$, 
$6.3 \times 10^{-3}$ and 
$2.3 \times 10^{-3}$, respectively. The distribution of these modes is shown in figure 6(b). The 
$\beta =0$ (spanwise uniform) mode is real and consists of a suction around the separation line and blowing around the backward stagnation point, similar to the 2-D optimal control to suppress vortex shedding (Mao et al. Reference Mao, Blackburn and Sherwin2015). The 
$\beta =1$ and 
$\beta =2$ (spanwise periodic) modes are concentrated around the separation line and reduce to trivial levels at the forward and backward stagnation points. This is slightly different from the spanwise periodic control to suppress instability (Guercio et al. Reference Guercio, Cossu and Pujals2014a), which peaks at 
$\theta ={\rm \pi} /2$ but has significant magnitude around the stagnation points. Each of the spanwise modes are symmetric or ‘in-phase’ with respect to the 
$x$ axis, which has been shown to be more effective in controlling the wake flow than the ‘out-of-phase’ profile (Kim & Choi Reference Kim and Choi2005; Guercio et al. Reference Guercio, Cossu and Pujals2014a).
4.2. Active control effects
The drag reduction effect of the active control is tested in DNS. The history of the drag coefficient, denoted as 
$\hat {C}_D$, at various values of the scale factor 
$S_a$ is presented in figure 7. We notice that the drag reduces monotonically with 
$|S_a|$ until 
$S_a=-7 \times 10^{-3}$, but the drag oscillation presents a more complicated dependence on 
$S_a$. At 
$S_a=-1 \times 10^{-3}$, the frequency of drag is almost unchanged, indicating a linear control effect. At 
$S_a=-3 \times 10^{-3}$, the oscillation at the frequency of uncontrolled vortex shedding is almost suppressed but a new oscillation with a lower frequency appears. At higher control magnitudes, this low-frequency oscillation is further suppressed, and at 
$S_a=-7 \times 10^{-3}$, the drag is almost time-independent. At even higher magnitude of the control, the drag does not further reduce and the oscillation reappears. Therefore, the following studies will be focused on 
$S_a=-7 \times 10^{-3}$. Shifting the initial condition of the uncontrolled flow in time or the spanwise direction and imposing the same control, similar drag reduction was observed, indicating that the reported control effect is insensitive to the initial condition to introduce the control.
Figure 7. Drag reduction at 
$Re=190$ and various control magnitudes. The control is scaled from sensitivity of drag to wall-normal active perturbations at 
$\omega =0$ (see figure 4a). The drag coefficient, denoted as 
$\hat {C}_D$, is averaged over the spanwise direction and will be used in all the following relevant plots.
Apart from 
$S_a$, the magnitude of the control can be also evaluated by the momentum coefficient, which is defined as the square integration of the control around the controlled boundary and the controlled time interval divided by the dynamic pressure of the uncontrolled flow. At 
$S_a=-7 \times 10^{-3}$, corresponding to maximum control magnitude 2 % of the free-stream velocity, the momentum coefficient is 
$10^{-4}$ and the drag reduction reaches 20 %. To achieve a similar drag reduction effect, the maximum control magnitude and momentum coefficient are above 8 % and 0.001, respectively, for both spanwise uniform suction and spanwise periodic blowing/suction reported in literature (Kim & Choi Reference Kim and Choi2005; Mao et al. Reference Mao, Blackburn and Sherwin2015).
The present control is suboptimized in both spanwise and azimuthal directions, while it has been shown that an optimization in the azimuthal direction only (with a prescribed spanwise wavenumber) lifts the control effect significantly in an effort to suppress instability in a cylinder flow at 
$Re\le 100$ (Guercio et al. Reference Guercio, Cossu and Pujals2014a). It is also worth noting that at 
$S_a=-7 \times 10^{-3}$ the force oscillation is very different from the unperturbed one, indicating that the perturbation magnitude is out of the linear regime and large enough to modify the stability characteristics of the ‘base’ flow, or more specifically stabilize the base flow, as will be detailed in § 4.4. The lift fluctuation (not shown here) was also suppressed by the control, similar to the drag fluctuation, further confirming the stabilizing of the flow and the suppression of flow unsteadiness.
It has been discussed in § 4.1 that the present control can be decomposed into the superposition of spanwise modes with 
$\beta =0$, 
$\beta =1$ and 
$\beta =2$, corresponding to boundary norm 
$2.2 \times 10^{-3}$, 
$6.3 \times 10^{-3}$ and 
$2.3 \times 10^{-3}$, respectively (see figure 6b). As both spanwise uniform (
$\beta =0$) and spanwise periodic (
$\beta \ne 0$) controls have been extensively studied, the control effect of each spanwise component of the present control is studied individually for comparison. From figure 8(a), all three individual modes have drag reduction effects, while the mode with 
$\beta =1$, which has a higher magnitude than the other two modes, reduces much more drag than the other two. The original control scaled from sensitivity, denoted as ‘full control’, has stronger drag reduction effects than each of its individual modes but its magnitude is also larger. To compare the drag reduction effectiveness of these modes, they are scaled to boundary norm 
$7 \times 10^{-3}$, the same as the full control. From figure 8(b), even at the same magnitude, the full control is still the most effective one, followed by 
$\beta =1$ and then 
$\beta =2$ and 
$\beta =0$. It is also tested that the 
$\beta =1$ (or 
$\beta =2$) complex Fourier mode has similar drag reduction effects as its real or imaginary part after scaling to the same boundary norm (not shown here). Therefore, over the parameters studied, the full control, in the form of a superposition of spanwise modes and localized in both spanwise and azimuthal directions, is more efficient than a single spanwise wave for drag reduction.
Figure 8. Drag reduction induced by individual spanwise modes of the control. In panel (a), the control is the spanwise-decomposed components of the scaled sensitivity of drag to wall-normal active perturbations at 
$Re=190$, 
$\omega =0$ and 
$S_a=-7 \times 10^{-3}$ (see figure 6b). In panel (b), all the individual modes are scaled to boundary norm 
$7 \times 10^{-3}$.
The active control is then applied to higher Reynolds numbers to test its robustness, as shown in figure 9. At 
$Re=1000$, the highest Reynolds number considered in this work, the uncontrolled drag coefficient is 1.0066, agreeing well with that reported in the literature by both numerical simulations and experiments (Wieselberger Reference Wieselberger1921; Henderson Reference Henderson1997). Similarly, with the 
$Re=190$ case, at 
$Re=300$ and 
$Re=1000$, significant drag reduction up to 20 % and the suppression of force oscillations can be observed. It is noted that at these higher Reynolds numbers, mode A becomes disordered but the control obtained at 
$Re=190$ is still effective, since it targets the 2-D shedding vortices. A similar result has been observed at 
$Re=300$ where mode B appears (Poncet et al. Reference Poncet, Hildebrand, Cottet and Koumoutsakos2008). These results suggest that in optimal control studies, when the sensitivity or gradient calculation at high 
$Re$ diverges owing to the chaotic dynamics of the flow (Wang & Gao Reference Wang and Gao2013), the sensitivity obtained at lower 
$Re$ can be a good approximation and yields significant control effects.
Figure 9. Controlled drag coefficients at various Reynolds numbers. The control is scaled from sensitivity of drag to wall-normal active perturbations at 
$\omega =0$, 
$S_a=-7 \times 10^{-3}$ and 
$Re=190$.
The benefit of the control can be evaluated by the ratio of the saved drag power and the actuation power to generate the control. For the present blowing/suction control, the dimensionless ideal actuation power can be defined as (Naito & Fukagata Reference Naito and Fukagata2014; Mao et al. Reference Mao, Blackburn and Sherwin2015)
\begin{equation} \dot{W}_B={-}\left\{ \hat{p}+\dfrac{1}{2}u_a^2+\dfrac{4}{Re}u_a, u_a \right\}, \end{equation}
where 
$\hat {p}$ is the controlled pressure and 
$u_a$ is the wall-normal control velocity (
$u_a \boldsymbol {n}=\boldsymbol {u}_a$). This idealized actuation power has been reported to be negative in a 2-D blowing/suction control of cylinder flow, where the actuators are distributed all around the cylinder surface with a maximum control velocity 0.4 (Naito & Fukagata Reference Naito and Fukagata2014). A negative power suggests that the control can be achieved by extracting power from the flow if all the actuators are ideally coordinated (the flow power extracted from one actuator can be applied to drive another actuator). Clearly such a coordination complicates the actuation system significantly and is impractical in real applications. Therefore, a maximum possible actuation power is defined to measure the control cost in a more practical manner:
\begin{equation} \dot{W}_{B\text{-}max}= \left\{ |\hat{p}|+\dfrac{1}{2}u_a^2+\dfrac{4}{Re}|u_a|, |u_a| \right\}. \end{equation}
For the control used above at 
$S_a=-7 \times 10^{-3}$, the ideal, maximum actuation power and drag coefficients are shown in table 2. The power savings ratio, defined as the ratio of control-reduced drag power and maximum input power, is over 20, indicating that the energy saved due to the control is over 20 times the energy to generate the control. This ratio is much higher than that of a 2-D suboptimal control, which is in the range 1.3–2.5 (Naito & Fukagata Reference Naito and Fukagata2014). It is worth noting that for surface rotation control, e.g. the perfect slip (free rotation of the surface driven by the flow), the actuator power is zero by definition (Shukla & Arakeri Reference Shukla and Arakeri2013).
Table 2. The uncontrolled drag coefficient (
$\hat {C}_D$ at 
$S_a=0$), controlled drag coefficient (
$\hat {C}_D$ at 
$S_a=-7 \times 10^{-3}$), ideal actuation power, maximum actuation power and power savings ratio at various Reynolds numbers. The control is scaled from the sensitivity of drag to wall-normal perturbations at 
$Re=190$ and 
$\omega =0$. The power savings ratio is defined as the ratio of the control-reduced drag and the maximum actuation power 
$\dot {W}_{B\text {-}max}$.
4.3. Drag reduction effect of the simplified control
To facilitate experimental implementations, the control obtained from (4.1) at 
${S_a=-7 \times 10^{-3}}$ is further simplified by retaining only the suction velocity above 
$0.005$, as shown in figure 10. This reduced control consists of two suctions localized in the spanwise and azimuthal directions and can be generated by only two actuators. Drag under this simplified control is shown in figure 11. Clearly this control produces a prominent drag reduction effect but weaker than that of the full control scaled from the sensitivity (see figure 9). At 
$Re=1000$, the drag reductions of the full control scaled from the sensitivity and the simplified control are 20 % and 15 %, respectively.
Figure 10. Distribution of the simplified control. Positive contour levels denote suction control.
Figure 11. Drag reduction under the simplified control in figure 10.
Then the actuation power is computed, similarly as performed for the full control. The two costs are 
$\dot {W}_B=5.2 \times 10^{-4}$ and 
$\dot {W}_{B\text {-}max}=4.1 \times 10^{-3}$ at 
$Re=1000$, corresponding to power savings ratio 17.8. Therefore, the simplified control still presents significant drag reduction effects at a cost below 6 % of the saved energy.
4.4. Mechanism of drag reduction
In this section, the flow fields are analysed to reveal the mechanism of drag reduction induced by the active control. Figure 12(a) illustrates the uncontrolled contour of instantaneous spanwise vorticity on the 
$z=0$ plane, while figures 12(b) and 12(c) show the controlled flow at 
$S_a=-5 \times 10^{-3}$ and 
$S_a=-7 \times 10^{-3}$, respectively. In the near wake, the averaged vortex sheets are not attached to the wall by the control as has been observed in optimal 2-D control, which delays boundary-layer separation, stabilizes the wake flow and subsequently suppresses vortex shedding (Mao Reference Mao2015a). Alternatively, the sheets are pushed away from each other and are attenuated, which suppresses and delays the vortex shedding. As the vortices are associated with local pressure minimum, marked by the black circles in the figure, the control pushes the low-pressure region further downstream away from the cylinder. Similar observations were obtained on other 
$z$ planes and they were not shown here for brevity.
Figure 12. Contours of the instantaneous spanwise vorticity on the 
$z=0$ plane. The black circles mark the points where pressure reaches local minimum. The control is scaled from the wall-normal sensitivity at 
$Re=190$ and 
$\omega =0$, and the magnitude of the control is 
$S_a=0$ (without control), 
$-5 \times 10^{-3}$ and 
$-7 \times 10^{-3}$ for (a,b,c), respectively.
The pressure distribution is further illustrated using the spanwise and temporally averaged flow, as shown in figure 13. It is clear that as the control becomes stronger from 
$S_a=0$ (without control) to 
$-7 \times 10^{-3}$, the pressure in the wake is significantly increased and the local pressure minimum is pushed away from the cylinder as observed in the instantaneous field in figure 12. The control has little impact on the boundary separation (the separation point is marked as back squares in figure 13) but enlarges the recirculation zone. Again these are opposite with 2-D optimal controls which suppress the boundary layer separation and the generation of recirculation bubbles, and imply a different control mechanism.
Figure 13. Streamlines and pressure contours in the spanwise and temporally averaged flow. The control is scaled from the wall-normal sensitivity at 
$Re=190$ and 
$\omega =0$, and the magnitude of the control is 
$S_a=0$ (without control), 
$-5 \times 10^{-3}$ and 
$-7 \times 10^{-3}$ for (a,b,c), respectively.
The streamlines in figure 13 illustrate a virtual surface covering the cylinder and the near wake recirculation zone, and the drag acting on this surface can be considered as the drag on the cylinder. When there is no control, this virtual surface is short and closes in the low-pressure region, and the flow around this surface is similar to flow around a bluff body at low Reynolds number without separation. When the control is imposed, this virtual surface is extended downstream and closes downstream of the region with local pressure minimum, and the flow around the surface resembles flow around a streamlined body.
The time and spanwise averaged pressure distribution on the surface of the cylinder is further plotted in figure 14. It is seen that the pressure is almost unmodified by the control around the fore part of the cylinder but is significantly increased downstream of the separation lines. This observation confirms that the control reduces drag by lifting the base pressure, a direct consequence of the suppression of wake oscillation. In this process, the attached boundary layer is not changed prominently, as has been seen in figures 12 and 13.
Figure 14. Spanwise and temporally averaged distribution of the pressure on the surface of the cylinder. In the controlled cases, the control magnitude is 
$S_a=-7 \times 10^{-3}$.
The 3-D structures are further examined to reveal the mechanism of vortex sheet attenuation. As shown in figure 15(a), in the near wake, the spanwise localized suction control distorts the vortex sheets towards the cylinder surface at around 
$z=2$, resulting in the localized vertical displacement of the other sheet. For example, the localized upward motion of the lower sheet owing to the suction effect on the surface pushes the upper sheet upwards, denoted as ‘
$y$-displacement’ in figures 15(a) and 15(b). The displaced part of the vortex sheet generates an induced velocity on the other parts, and causes the streamwise stretch of the overall vortex sheet, marked as ‘
$x$-stretch’ in figure 15(c). Clearly the ‘
$y$-displacement’ and the ‘
$x$-stretch’ reduce the interaction of the vortex sheets by pushing them away from each other, and attenuate their strength in the spanwise averaged sense (also see figure 12). In this process, the flow becomes more inhomogeneous in the spanwise direction, which can be observed by comparing the vortex lines shown in figures 3 and 15(a).
Figure 15. Iso-surfaces of spanwise vorticity 0.5 (red) and 
$-0.5$ (blue) at 
$S_a=-7 \times 10^{-3}$ and 
$Re=190$. Panels (a,b,c) are 3-D, side (
$x$–
$y$ plane) and top (
$x$–
$z$ plane) views of the same surface, respectively. The flow field is repeated in 
$z$ to highlight the spanwise variation. The thick solid lines in panel (a) represent the vortex lines associated with the lower vortex sheet.
The development of the control induced flow is schematically plotted. Figure 16(a) shows the flow patten at 
$z=0$ or the uncontrolled vortex shedding scenario. Figure 16(b) illustrates the effects of spanwise localized (around 
$z=2$) suction on the vortex shedding, which was displaced in the vertical direction to be away from the 
$x$ axis. Then figure 16(c) shows that the displaced part of the vortex sheet generates an induced velocity in the 
$x$ direction to the other parts of the sheet. This induced velocity stretches the vortex sheet in the 
$x$ direction and attenuates it significantly as illustrated in figure 16(d). It is noted that this control mechanism is robust and is observed at other control levels and Reynolds numbers.
Figure 16. Schematics of the drag reduction mechanism. Panel (a) is the vortex shedding at 
$z=0$ showing the uncontrolled scenario; panel (b) shows the suction effect of the control and its induced 
$y$-displacement; panel (c) illustrates the displacement induced velocity in the streamwise direction; panel (d) presents the stretch in streamwise direction owing to the induced velocity. Lengths of the arrows in panel (c) represent the magnitude of the induced velocity. Blue and red denote negative and positive vorticity, respectively.
Therefore, the drag reduction effect observed in the present work stems from both a suppression of vortical structure interactions via the 
$y$-displacement, and the attenuation of the vortex structures via the 
$x$-stretch. Similar stretch of the wake structure and its induced spanwise mismatch has been observed before in drag reduction and instability investigations using spanwise periodic forcing (Kim & Choi Reference Kim and Choi2005; Hwang et al. Reference Hwang, Kim and Choi2013). This attenuation of wake structures pushes the vortex shedding associated with local pressure minimum away from the cylinder to generate a drag reduction effect (Poncet et al. Reference Poncet, Hildebrand, Cottet and Koumoutsakos2008). In the controlled flow, the separation is not postponed and the recirculation zone is enlarged. This control mechanism is clearly different from those associated with suppressing separation and attenuating the recirculation by means of either blowing/suction or surface slip mainly in two dimensions (e.g. Min & Choi Reference Min and Choi1999; Delaunay & Kaiktsis Reference Delaunay and Kaiktsis2001; Legendre, Laura & Magnaudet Reference Legendre, Laura and Magnaudet2009; Shukla & Arakeri Reference Shukla and Arakeri2013; Mao Reference Mao2015a) and with some 3-D exceptions (e.g. Shtendel & Seifert Reference Shtendel and Seifert2014).
5. Passive drag reduction
In this section, the sensitivity to passive perturbations and the corresponding passive control effect on drag will be examined. Firstly, the sensitivity of drag with respect to passive perturbations is computed, as show in figure 17(a). This sensitivity consists of spanwise uniform and periodic components both concentrated around the separation lines, similar to the sensitivity to active perturbations shown in figure 4.
Figure 17. Panel (a) is the contour of the sensitivity to wall-normal passive perturbations normalized by its boundary norm at 
$Re=190$. The dash-dotted lines represent the mean separation lines. Panel (b) is the surface of the cylinder with roughness scaled from the sensitivity at 
$S_r=-0.04$.
Then a passive control can be obtained by scaling this sensitivity:
\begin{equation} \boldsymbol{B}_r=\boldsymbol{B}+ S_r \frac{\boldsymbol{\nabla}_{r}\bar{f}_{\boldsymbol{K}}}{ \|\boldsymbol{\nabla}_{r}\bar{f}_{\boldsymbol{K}}\|}, \end{equation}
where 
$S_r$ is a scale factor and its absolute value denotes the boundary norm of the passive perturbation or surface roughness. The surface of the cylinder with passive control at 
$S_r=-0.04$ is illustrated in figure 17(b).
The passive control with a single spanwise wavenumber has been observed to suppress vortex shedding effectively; surface waviness of maximum height 1.5 % of the cylinder diameter is sufficient to stabilize the wake oscillation at 
$Re=100$ (Tammisola Reference Tammisola2017). We note that unlike the spanwise uniform or periodic passive control, the present one presents a ‘concentrated roughness’, where the surface deformation only occurs over a small portion of the surface in both azimuthal and spanwise directions. Therefore, the implementation of this control only requires a localized modification of the geometry.
The drag history at various 
$S_r$ is shown in figure 18. Data was collected from 
$t=1000$ to exclude the initial transient period. The drag history at increasing magnitudes of the passive control is similar to the active control case. From 
$S_r=-0.01$ to 
$S_r=-0.04$, both the drag and the drag oscillation reduce, indicating a stabilization of the flow field. At 
$S_r=-0.04$, the drag reduction reaches 13 %. At even larger magnitude of the passive control 
$S_r=-0.06$, the drag is still reducing but the oscillation becomes prominent. As 
$S_r=-0.06$ corresponds to a significant change of the shape of the cylinder (maximum height of roughness over 0.12), it will not be analysed in the following.
Figure 18. Passively controlled drag coefficient at 
$Re=190$. The control is scaled from the passive sensitivity of drag to wall-normal roughness. The surface with passive control at 
$S_r=-0.04$ has been illustrated in figure 17(b).
Similarly, as in the active control studies, this passive control or surface deformation is also tested at higher Reynolds numbers. As shown in figure 19, the drag reduction is still significant up to 
$Re=1000$, the highest Reynolds number considered. The mechanism of drag reduction is found to be similar to that observed for the active control consisting of a 
$y$-displacement and an 
$x$-stretch, and is not repeated here.
Figure 19. Drag coefficients under passive control (see figure 17(b) for the deformed geometry) at various Reynolds numbers. The control is scaled from the sensitivity of drag to wall-normal roughness at 
$S_r=-0.04$ and 
$Re=190$.
6. Conclusion
The well-established algorithm to compute sensitivities of modal and non-modal instabilities to perturbations is extended to calculate the sensitivity of forces to surface control forcing in flow past a cylinder. Both active and passive controls, regarded as boundary velocity perturbations and surface roughness, respectively, are considered. The sensitivity with respect to active perturbations had been established and applied to 2-D uncontrolled flow (Mao Reference Mao2015b) while the passive counterpart is developed in the present work.
The algorithms are then applied to flow past a circular cylinder at Reynolds number 
$Re=190$, where the flow field is three-dimensional. Instead of assigning a prescribed spanwise wavenumber to the perturbations as documented in the literature (Hwang et al. Reference Hwang, Kim and Choi2013; Guercio et al. Reference Guercio, Cossu and Pujals2014a,Reference Guercio, Cossu and Pujalsb; Tammisola et al. Reference Tammisola, Giannetti, Citro and Juniper2014), the surface perturbation, as well as the sensitivity, is represented as a superposition of multiple spanwise Fourier modes.
The sensitivity is further scaled as a control to reduce drag. The active control is localized in both spanwise and azimuthal directions, and can be decomposed as waves with spanwise wavenumber 
$\beta =0$, 
$\beta =1$ and 
$\beta =2$. It can be expected that at even higher Reynolds numbers, this control will be more localized and therefore can be generated by a limited number of blowing/suction actuators in experiments. At maximum control magnitude 2 % of the free-stream velocity (or momentum coefficient 
$10^{-4}$), the drag is reduced by 20 % and the temporal oscillation of drag is almost suppressed. This control is more effective than each of its spanwise modes with or without scaling the modes to the same magnitude as the control. To achieve a similar drag reduction effect, a maximum control magnitude above 8 % (or momentum coefficient above 
$10^{-3}$) is required for either the uniform suction or periodic blowing/suction reported in literature. The scaled sensitivity obtained at 
$Re=190$ is then applied to flow at higher 
$Re$. The power savings ratio, defined as the ratio of control-reduced drag power and maximum input power, is above 
$20$, in the range 
$190\le Re\le 1000$. The control is further simplified to spanwise localized suctions around the separation lines to facilitate experimental implementations. This simplified control still produces significant drag reduction effects but weaker than that of the full control. The passive control is also localized in azimuthal and spanwise directions, and therefore its implementation only requires a localized change of the surface geometry. Significant drag reduction is also observed in the passively controlled cylinder flow at 
$Re\le 1000$ and the mechanism is similar to that observed in active control.
Examining the controlled flow field, it is noticed that the localized suction control induces a 
$y$-displacement and an 
$x$-stretch to the vorticity sheets and attenuates their strength. The attenuation of vortex structures has been observed to stabilize the flow and reduce drag in spanwise periodic control (Kim & Choi Reference Kim and Choi2005; Hwang et al. Reference Hwang, Kim and Choi2013; Guercio et al. Reference Guercio, Cossu and Pujals2014b; Tammisola Reference Tammisola2017), while the present study indicates that a localized control, which facilitates physical implementations, is more effective to generate this effect. As a consequence of the displacement and stretch, the vortex shedding is delayed and the local pressure minimum is pushed away from the cylinder, which induces a drag reduction effect. Different from the 2-D control which suppresses boundary-layer separation and the generation of recirculation bubbles (Min & Choi Reference Min and Choi1999; Delaunay & Kaiktsis Reference Delaunay and Kaiktsis2001; Legendre et al. Reference Legendre, Laura and Magnaudet2009; Shukla & Arakeri Reference Shukla and Arakeri2013; Mao Reference Mao2015a), the present 3-D control enlarges the bubble and produces a virtual surface to convert the cylinder flow to a streamlined body flow.
Funding
The authors would like to acknowledge financial support under Engineering and Physical Sciences Research Council (EPSRC) grant EP/M025039/2. Calculations were performed on HPC-Midlands and ARCHER funded by the EPSRC (grant number EP/K000063/1 and EP/R029326/1, respectively).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Adjoint analyses
Considering integration by parts and the divergence theorem (Barkley et al. Reference Barkley, Blackburn and Sherwin2008), there is
\begin{align} -\langle \boldsymbol{u}^*, \partial_t \boldsymbol{u} - {\mathcal{L}}(\boldsymbol{u})\rangle = \langle \boldsymbol{u}, \partial_t \boldsymbol{u}^* + {\mathcal{L}}^*(\boldsymbol{u}^*)\rangle - \tau ^{{-}1} (\boldsymbol{u}_{\tau}, \boldsymbol{u}^*_{\tau}) +\tau ^{{-}1} (\boldsymbol{u}_0, \boldsymbol{u}^*_0) \nonumber\\ +\tau ^{{-}1} \int_0^\tau \int_{\partial\varOmega}\boldsymbol{n}\boldsymbol{\cdot}[-\boldsymbol{U}(\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{u}^*) + \boldsymbol{u}p^* - \boldsymbol{u}^*p + Re^{{-}1} (\boldsymbol{\nabla} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{u}^*-\boldsymbol{\nabla} \boldsymbol{u}^*\boldsymbol{\cdot}\boldsymbol{u})] \, \mathrm{d} \partial \varOmega \,\mathrm{d} t, \end{align}
where 
$\partial \varOmega$ represents all the boundaries of the computational domain; superscript 
$*$ denotes the adjoint variables; 
$\boldsymbol {u}^*_{\tau }$ and 
$\boldsymbol {u}_{\tau }$ are the adjoint velocity and perturbation velocity at 
$t=\tau$, respectively. Here 
${\mathcal {L}}$ is an operator and 
$\partial _t \boldsymbol {u}-{\mathcal {L}}(\boldsymbol {u}) =0$ represents the linearized NS equation (2.2), 
${\mathcal {L}}^*$ is the adjoint operator of 
${\mathcal {L}}$ and correspondingly 
$\partial _t \boldsymbol {u}^* + {\mathcal {L}}^*(\boldsymbol {u}^*)=0$ is the adjoint of the linearized NS equation, which can be expanded as
\begin{equation} \partial_t \boldsymbol{u}^*+ \boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^*- \boldsymbol{\nabla} \boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{u}^*- \boldsymbol{\nabla} p^*+ {\textit{Re}}^{{-}1} \nabla^2 \boldsymbol{u}^*=0 \quad \mathrm{and} \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}^*=0. \end{equation}
On the inflow and far-field boundaries, zero-Dirichlet and high-order-Neumann conditions are used for adjoint velocity and pressure, respectively (Karniadakis, Israeli & Orszag Reference Karniadakis, Israeli and Orszag1991); on the outflow, a mixed velocity boundary condition 
${{\textit {Re}} ^{-1}\partial _{\boldsymbol {n}}\boldsymbol {u}^*+ \boldsymbol {n} \boldsymbol {\cdot } (\boldsymbol {U} +\boldsymbol {u}) \boldsymbol {u}^*=0}$ and zero-Dirichlet pressure condition are implemented (Mao et al. Reference Mao, Blackburn and Sherwin2013); on the perturbed boundary, a Dirichlet velocity and computed Neumann pressure conditions are used. Considering the signs of the viscous term and of the time derivative term, this adjoint equation should be initialized at 
$t=\tau$ and integrated backwards to 
$t=0$.
Since 
$\boldsymbol {u}$ satisfies the linearized NS equation and its initial condition is zero, 
${\langle \boldsymbol {u}^*, \partial _t \boldsymbol {u} - {\mathcal {L}}(\boldsymbol {u})\rangle =0}$ and 
${(\boldsymbol {u}_0, \boldsymbol {u}^*_0)=0}$. Then if the adjoint variables are solutions of the adjoint equation initialized by 
$\boldsymbol {u}_{\tau }^*=0$, all the terms in (A1) are zero except the last one. The choice of boundary conditions for the adjoint variables ensures that the integration over the inflow, outflow and far-field boundaries in this last term is zero. Therefore, one obtains
\begin{equation} \{ p \boldsymbol{n} - Re^{{-}1} \partial_{\boldsymbol{n}}\boldsymbol{u}, \boldsymbol{u}^* \} =\{ p^* \boldsymbol{n} - Re^{{-}1}\partial_{\boldsymbol{n}}\boldsymbol{u}^*,\boldsymbol{u} \}. \end{equation}
