2.1. Modelling the unsteadiness

To simplify notation, we introduce scalar products

(2.1a−c ) $$\begin{eqnarray}(\boldsymbol{a},\boldsymbol{b})=\int _{{\it\Omega}}\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{b}\,\text{d}V,\quad \langle \boldsymbol{c},\boldsymbol{d}\rangle ={\it\tau}^{-1}\int _{0}^{{\it\tau}}(\boldsymbol{c},\boldsymbol{d})\,\text{d}t,\quad [\boldsymbol{e},\boldsymbol{f}]=\int _{\partial {\it\Omega}_{c}}\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{f}\,\text{d}S,\end{eqnarray}$$

where ${\it\Omega}$ represents the spatial domain, $\partial {\it\Omega}_{c}$ denotes the control boundary, referring to the boundary segment where the control is imposed, $t$ is time, ${\it\tau}$ is a final time, $\boldsymbol{a}$ and $\boldsymbol{b}$ are vector fields defined on the spatial domain ${\it\Omega}$ , $\boldsymbol{c}$ and $\boldsymbol{d}$ are defined on the spatial domain ${\it\Omega}$ and time domain $[0,{\it\tau}]$ , and $\boldsymbol{e}$ and $\boldsymbol{f}$ are defined on the control boundary $\partial {\it\Omega}_{c}$ .

In an initial test of the control effects of boundary perturbations, it is observed that the optimal control is concentrated on the wall-normal component. Therefore, in the following, we consider Dirichlet-type boundary-normal velocity perturbations as the control, denoted as $\boldsymbol{u}_{c}(\boldsymbol{x},t)$ , where $\boldsymbol{x}$ represents the coordinate of the control boundary. To reduce the dimension of $\boldsymbol{u}_{c}(\boldsymbol{x},t)$ after temporal–spatial discretization, we decompose the temporal and spatial dependence as

(2.2) $$\begin{eqnarray}\boldsymbol{u}_{c}(\boldsymbol{x},t)=u_{n}(\boldsymbol{x})f(t,{\it\omega})\boldsymbol{n},\end{eqnarray}$$

where $u_{n}(\boldsymbol{x})$ is the spatial dependence of the perturbation on the wall-normal direction to be optimized, $\boldsymbol{n}$ is the outward-pointing unit normal vector on the boundary and $f(t,{\it\omega})$ is a prescribed temporal-dependence function. In this work, we adopt a temporal function

(2.3) $$\begin{eqnarray}f(t,{\it\omega})=[1-\exp (-t^{2})]\cos ({\it\omega}t),\end{eqnarray}$$

where the first term on the right ensures $\boldsymbol{u}_{c}(\boldsymbol{x},0)=0$ so as to avoid temporal discontinuity in the numerical simulation at $t=0$ . When the final time ${\it\tau}$ is large enough, ${\it\omega}$ becomes the frequency of the control perturbation.

The flow field can be decomposed as the sum of an uncontrolled base flow and a perturbation flow induced by the control, i.e.  $(\hat{\boldsymbol{u}},\hat{p})=(\boldsymbol{U},P)+(\boldsymbol{u},p)$ , where $\hat{\boldsymbol{u}}$ , $\boldsymbol{U}$ and $\boldsymbol{u}$ are respectively the total velocity, uncontrolled velocity and perturbation velocity, while $\hat{p}$ , $P$ and $p$ are the total pressure, uncontrolled pressure and perturbation pressure. Then, we substitute this decomposition into the NS equation to reach

(2.4) $$\begin{eqnarray}\partial _{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}+\boldsymbol{{\rm\nabla}}p-\mathit{Re}^{-1}\boldsymbol{{\rm\nabla}}^{2}\boldsymbol{u}=0,\quad \text{with}~\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0.\end{eqnarray}$$

In the following, this equation is written in a compact form as

(2.5) $$\begin{eqnarray}\partial _{t}\boldsymbol{u}-L(\boldsymbol{u})+N(\boldsymbol{u})=0,\end{eqnarray}$$

where $N(\boldsymbol{u})=\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}$ denotes the nonlinear convection term and $L(\boldsymbol{u})$ represents the pressure and viscous terms as well as the divergence-free condition, which are linear. In the following this equation will be referred to as the decomposed NS equation.

Since the control perturbation is imposed at $t=0$ , the initial condition of perturbation $\boldsymbol{u}$ is set to zero. This initial condition is compatible with the control perturbation defined in (2.2).

On the inflow and far-field boundaries, zero Dirichlet and computed Neumann conditions are adopted for $\boldsymbol{u}$ and $p$ (Blackburn & Sherwin Reference Blackburn and Sherwin2004) respectively; on the outflow boundary, $\boldsymbol{{\rm\nabla}}_{\boldsymbol{n}}\boldsymbol{u}=0$ and $p=0$ are imposed; on the control boundary, the Dirichlet boundary condition and computed Neumann condition are implemented for $\boldsymbol{u}$ and $p$ respectively.

For convenience, we define an operator $\mathscr{N}$ corresponding to the action of evolving $u_{n}(\boldsymbol{x})$ to $\boldsymbol{u}(t)$ by integrating the decomposed NS equation (2.5), i.e. 

(2.6) $$\begin{eqnarray}\boldsymbol{u}(t)=\mathscr{N}(u_{n}).\end{eqnarray}$$

The unsteadiness of the flow field over the time horizon $[0,{\it\tau}]$ can be modelled as the fluctuation energy

(2.7) $$\begin{eqnarray}A(\boldsymbol{u})=\langle W\hat{\boldsymbol{u}},W\hat{\boldsymbol{u}}\rangle -(W\hat{\boldsymbol{u}}_{a},W\hat{\boldsymbol{u}}_{a}),\end{eqnarray}$$

where $W$ is a non-negative spatial weight used to filter the unsteadiness in the regions not of interest when we are only concerned about the unsteadiness in a localized region of the computational domain. When the unsteadiness across the whole domain is taken into account uniformly, $W(\boldsymbol{x})=1$ . Here, $\hat{\boldsymbol{u}}_{a}$ is the averaged velocity, calculated as

(2.8a,b ) $$\begin{eqnarray}\hat{\boldsymbol{u}}_{a}=\boldsymbol{U}_{a}+\boldsymbol{u}_{a},\quad \text{where}~\boldsymbol{U}_{a}={\it\tau}^{-1}\int _{0}^{{\it\tau}}\boldsymbol{U}\,\text{d}t\quad \text{and}\quad \boldsymbol{u}_{a}={\it\tau}^{-1}\int _{0}^{{\it\tau}}\boldsymbol{u}\,\text{d}t\end{eqnarray}$$

represent the averaged velocities of base flow and perturbation flow respectively. For flow around a cylinder at the relatively low Reynolds number considered in this work, this unsteadiness quantifies the magnitude of vortex shedding, while for a fully turbulent flow at higher Reynolds number, it describes the turbulence intensity.

2.2. Lagrangian functional

Similarly to the methodology established by Mao et al. (Reference Mao, Blackburn and Sherwin2012) to calculate the optimal boundary perturbation that results in maximum energy in the domain over a finite time horizon, a Lagrangian functional can be defined as

(2.9) $$\begin{eqnarray}\mathscr{L}=A-\langle \boldsymbol{u}^{\ast },\partial _{t}\boldsymbol{u}-L(\boldsymbol{u})+N(\boldsymbol{u})\rangle +{\it\lambda}(E-[u_{n},u_{n}]),\end{eqnarray}$$

where $\boldsymbol{u}^{\ast }$ represents the adjoint velocity, ${\it\lambda}$ is a Lagrange multiplier and $E$ is a prescribed control cost which is equal to the square of the wall transpiration velocity integrated over the controlled surface. Hence, for the circular cylinder, an equivalent root-mean-square (r.m.s.) value of the transpiration velocity may be computed as $(E/{\rm\pi})^{1/2}$ . The first term following the equality in (2.9) represents the unsteadiness or the fluctuation energy to be minimized, the second term applies the constraint of the governing equation, which is the decomposed NS equation (2.5), and the third term is the constraint on the control cost or the magnitude of the control perturbation. Various options for the constraint on the magnitude of perturbations have been discussed by Foures, Caulfield & Schmid (Reference Foures, Caulfield and Schmid2012), while the one we adopt reflects the momentum coefficient.

The adjoint variables are calculated from an adjoint equation. To facilitate the derivation of the adjoint equation, we reformulate the Lagrangian functional. Using integration by parts and the divergence theorem (Barkley, Blackburn & Sherwin Reference Barkley, Blackburn and Sherwin2008), we have

(2.10) $$\begin{eqnarray}\displaystyle -\langle \boldsymbol{u}^{\ast },\partial _{t}\boldsymbol{u}-L(\boldsymbol{u})\rangle & = & \displaystyle \langle \boldsymbol{u},\partial _{t}\boldsymbol{u}^{\ast }+L^{\ast }(\boldsymbol{u}^{\ast })\rangle -{\it\tau}^{-1}(\boldsymbol{u}_{{\it\tau}},\boldsymbol{u}_{{\it\tau}}^{\ast })\nonumber\\ \displaystyle & & \displaystyle +\,{\it\tau}^{-1}\int _{0}^{{\it\tau}}\int _{\partial {\it\Omega}}\boldsymbol{n}\boldsymbol{\cdot }[-\boldsymbol{U}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\ast })+\boldsymbol{u}p^{\ast }-\boldsymbol{u}^{\ast }p\nonumber\\ \displaystyle & & \displaystyle +\,\mathit{Re}^{-1}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }-\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\boldsymbol{u})]\,\text{d}S\,\text{d}t,\end{eqnarray}$$

where $\partial {\it\Omega}$ represents all the boundaries of the computational domain; $\boldsymbol{u}_{{\it\tau}}^{\ast }$ and $\boldsymbol{u}_{{\it\tau}}$ are the adjoint and perturbation velocity vector at $t={\it\tau}$ respectively; $L^{\ast }(\boldsymbol{u}^{\ast })=\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\ast }-\boldsymbol{{\rm\nabla}}\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }-\boldsymbol{{\rm\nabla}}p^{\ast }+\mathit{Re}^{-1}\boldsymbol{{\rm\nabla}}^{2}\boldsymbol{u}^{\ast }$ and $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }=0$ . We note that $\partial _{t}\boldsymbol{u}^{\ast }+L^{\ast }(\boldsymbol{u}^{\ast })=0$ is the adjoint equation extensively used in investigations of receptivity and non-normality (Hill Reference Hill1995; Giannetti & Luchini Reference Giannetti and Luchini2007; Barkley et al. Reference Barkley, Blackburn and Sherwin2008). In the derivation of (2.10), the zero initial condition of $\boldsymbol{u}$ is implied. Substituting (2.10) into (2.9), we have a reformulated Lagrangian functional

(2.11) $$\begin{eqnarray}\displaystyle \mathscr{L} & = & \displaystyle A-\langle \boldsymbol{u}^{\ast },N(\boldsymbol{u})\rangle +{\it\lambda}(E-[u_{n},u_{n}])+\langle \boldsymbol{u},\partial _{t}\boldsymbol{u}^{\ast }+L^{\ast }(\boldsymbol{u}^{\ast })\rangle -{\it\tau}^{-1}(\boldsymbol{u}_{{\it\tau}},\boldsymbol{u}_{{\it\tau}}^{\ast })\nonumber\\ \displaystyle & & \displaystyle +\,{\it\tau}^{-1}\int _{0}^{{\it\tau}}\int _{\partial {\it\Omega}}\boldsymbol{n}\boldsymbol{\cdot }[-\boldsymbol{U}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\ast })+\boldsymbol{u}p^{\ast }-\boldsymbol{u}^{\ast }p\nonumber\\ \displaystyle & & \displaystyle +\,\mathit{Re}^{-1}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }-\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\boldsymbol{u})]\,\text{d}S\,\text{d}t.\end{eqnarray}$$

As stated above, the adjoint variables are calculated from an adjoint equation, which can be obtained by setting the first variation of the Lagrangian with respect to $\boldsymbol{u}$ to zero. On the right-hand side of (2.11), the variation of the first term with respect to the velocity vector $\boldsymbol{u}$ is

(2.12) $$\begin{eqnarray}{\it\delta}A({\it\delta}\boldsymbol{u})=\langle {\it\delta}\boldsymbol{u},2W^{2}(\boldsymbol{U}+\boldsymbol{u}-\boldsymbol{U}_{a}-\boldsymbol{u}_{a})\rangle ,\end{eqnarray}$$

where we have used the relation $(\boldsymbol{u}_{a},\boldsymbol{u}_{a})=\langle \boldsymbol{u}_{a},\boldsymbol{u}\rangle$ and similar manipulations. Next, we consider the variation of the second term with respect to $\boldsymbol{u}$ , i.e. 

(2.13) $$\begin{eqnarray}{\it\delta}\langle \boldsymbol{u}^{\ast },N(\boldsymbol{u})\rangle ({\it\delta}\boldsymbol{u})=\langle {\it\delta}\boldsymbol{u},\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }-\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\ast }\rangle +{\it\tau}^{-1}\int _{0}^{{\it\tau}}\int _{\partial {\it\Omega}}({\it\delta}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\ast })(\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{u})\,\text{d}S\,\text{d}t.\end{eqnarray}$$

(For details of this derivation, refer to appendix A.)

Combining (2.11)–(2.13), we have the variation of the Lagrangian with respect to  $\boldsymbol{u}$ ,

(2.14) $$\begin{eqnarray}\displaystyle {\it\delta}\mathscr{L}({\it\delta}\boldsymbol{u}) & = & \displaystyle \langle {\it\delta}\boldsymbol{u},\partial _{t}\boldsymbol{u}^{\ast }+L^{\ast }(\boldsymbol{u}^{\ast })-\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\ast }+2W^{2}(\boldsymbol{U}+\boldsymbol{u}-\boldsymbol{U}_{a}-\boldsymbol{u}_{a})\rangle \nonumber\\ \displaystyle & & \displaystyle -\,{\it\tau}^{-1}({\it\delta}\boldsymbol{u}_{{\it\tau}},\boldsymbol{u}_{{\it\tau}}^{\ast })+{\it\tau}^{-1}\int _{0}^{{\it\tau}}\int _{\partial {\it\Omega}}\boldsymbol{n}\boldsymbol{\cdot }[-(\boldsymbol{U}+\boldsymbol{u})({\it\delta}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\ast })+{\it\delta}\boldsymbol{u}p^{\ast }-\boldsymbol{u}^{\ast }{\it\delta}p({\it\delta}\boldsymbol{u})\nonumber\\ \displaystyle & & \displaystyle +\,\mathit{Re}^{-1}(\boldsymbol{{\rm\nabla}}{\it\delta}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }-\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\ast }\boldsymbol{\cdot }{\it\delta}\boldsymbol{u})]\,\text{d}S\,\text{d}t,\end{eqnarray}$$

where ${\it\delta}p({\it\delta}\boldsymbol{u})$ is the variation of $p$ induced by the variation of $\boldsymbol{u}$ . In this derivation, $\boldsymbol{u}$ , $\boldsymbol{u}^{\ast }$ and $u_{n}$ are considered as independent variables. Setting this variation to zero, we obtain the adjoint equation, its initial condition and boundary conditions.

From the first term on the right-hand side of (2.14) we have the adjoint equation:

(2.15) $$\begin{eqnarray}\partial _{t}\boldsymbol{u}^{\ast }+L^{\ast }(\boldsymbol{u}^{\ast })-\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\ast }+2W^{2}(\boldsymbol{U}+\boldsymbol{u}-\boldsymbol{U}_{a}-\boldsymbol{u}_{a})=0.\end{eqnarray}$$

Inspecting the sign of the time derivative term and diffusion term, we see that this equation should be integrated backwards from $t={\it\tau}$ to $t=0$ . The initial condition $\boldsymbol{u}_{{\it\tau}}^{\ast }$ can be obtained from the second term on the right-hand side of (2.14), that is $\boldsymbol{u}_{{\it\tau}}^{\ast }=0$ .

The boundary condition of this adjoint equation can be obtained from the last term of (2.14). The combination of adjoint velocity and pressure boundary conditions has to ensure that the integration in the last term of (2.14) over each boundary is zero. On the inflow boundary, control boundary and far-field boundary, the velocity is fixed and therefore ${\it\delta}\boldsymbol{u}=0$ . On these three boundaries, zero Dirichlet and computed Neumann conditions are used for adjoint variables $\boldsymbol{u}^{\ast }$ and $p^{\ast }$ respectively. On the outflow boundary, the boundary conditions for adjoint variables are $\mathit{Re}^{-1}\boldsymbol{{\rm\nabla}}_{\boldsymbol{n}}\boldsymbol{u}^{\ast }+\boldsymbol{n}\boldsymbol{\cdot }(\boldsymbol{U}+\boldsymbol{u})\boldsymbol{u}^{\ast }=0$ and $p^{\ast }=0$ , where $\boldsymbol{{\rm\nabla}}_{\boldsymbol{n}\boldsymbol{u}}^{\ast }=\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\ast }$ (Mao, Blackburn & Sherwin Reference Mao, Blackburn and Sherwin2013). The choices of these boundary conditions ensure that the integration over the boundaries in the last term in (2.14) is zero. It is worth noting that the adjoint velocity outflow boundary condition is dependent on the velocity vector $\boldsymbol{u}$ , and so this boundary condition has to be updated (to the latest $\boldsymbol{u}$ ) at every time step when solving the adjoint equation.

From the last term of (2.11), we have the variation of the Lagrangian with respect to the control perturbation $u_{n}$ ,

(2.16) $$\begin{eqnarray}{\it\delta}\mathscr{L}({\it\delta}u_{n})=[g-2{\it\lambda}u_{n},{\it\delta}u_{n}],\end{eqnarray}$$

where

(2.17) $$\begin{eqnarray}g={\it\tau}^{-1}\boldsymbol{n}\boldsymbol{\cdot }\int _{0}^{{\it\tau}}(p^{\ast }\boldsymbol{n}-\mathit{Re}^{-1}\boldsymbol{{\rm\nabla}}_{\boldsymbol{ n}}\boldsymbol{u}^{\ast })f(t,-{\it\omega})\,\text{d}t.\end{eqnarray}$$

Here, $f(t,-{\it\omega})$ is the adjoint of $f(t,{\it\omega})$ , satisfying $[\,f(t,{\it\omega})\boldsymbol{e},\boldsymbol{f}]=[\boldsymbol{e},f(t,-{\it\omega})\boldsymbol{f}]$ , with $\boldsymbol{e}$ and $\boldsymbol{f}$ representing vectors defined on the control boundary.

From the definition of the gradient of the Lagrangian associated with the Gâteaux differential, the gradient of $\mathscr{L}$ with respect to $u_{n}$ is

(2.18) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}_{\boldsymbol{u}_{\boldsymbol{n}}}\mathscr{L}=g-2{\it\lambda}u_{n}.\end{eqnarray}$$

In this formulation, both $g$ and $u_{n}$ are functions defined on the control boundary.

2.3. Nonlinear optimal control

To compute an optimal solution of $u_{n}$ , an iterative optimization algorithm is required. For the gradient-based optimization methods, in each iteration, the gradient of the Lagrangian with respect to the control has to be evaluated. From the expression of the gradient in (2.18), we see that the adjoint equation must be solved to calculate $g$ and then the gradient. From the adjoint equation and its boundary conditions, we note that the full development history of $\boldsymbol{u}$ is required to solve the adjoint variables. Solving the adjoint equation is computationally as intense as solving the decomposed NS equation, but requires much larger memory. To reduce this high requirement on memory, a checkpointing scheme has been proposed by Griewank (Reference Griewank1992). In that scheme, the memory requirement can be significantly reduced but the governing equation (2.5) and the adjoint equation (2.15) have to be restarted at every checkpoint. In the numerical method we presently implement to solve these equations, large-scale matrices need to be generated and installed every time the equation solvers are restarted, and therefore such a checkpointing scheme is inefficient. In the current work, the development history of $\boldsymbol{u}$ at every time step is written to file when integrating the decomposed NS equation (2.5) and read into memory at every step with the reversed order when integrating the adjoint equation (2.15).

In the optimization algorithm, a steepest-descent method is used and therefore the search direction is the gradient itself. The control perturbation and perturbation velocity are updated along this direction from iterative step $k$ to $k+1$ as

(2.19) $$\begin{eqnarray}u_{n}^{k+1}=u_{n}^{k}+{\it\alpha}\boldsymbol{{\rm\nabla}}_{\boldsymbol{ u}_{\boldsymbol{n}}}\mathscr{L}^{k}=u_{n}^{k}+{\it\alpha}g^{k}-2{\it\alpha}{\it\lambda}u_{n}^{k},\end{eqnarray}$$

where ${\it\alpha}$ is a step length and ${\it\lambda}$ can be obtained by imposing the constraint of control cost on $u_{n}^{k+1}$ .

Correspondingly, the perturbation velocity at step $k+1$ can be obtained as

(2.20) $$\begin{eqnarray}\boldsymbol{u}^{k+1}(t)=\mathscr{N}(u_{n}^{k}+{\it\alpha}g^{k}-2{\it\alpha}{\it\lambda}u_{n}^{k}).\end{eqnarray}$$

If the step length ${\it\alpha}$ is small enough, ${\it\alpha}g^{k}$ and $-2{\it\alpha}{\it\lambda}u_{n}^{k}$ are small compared with $u_{n}^{k}$ . Therefore, we have

(2.21) $$\begin{eqnarray}\boldsymbol{u}_{l}^{k+1}=\boldsymbol{u}^{k}+{\it\alpha}\mathscr{M}_{\boldsymbol{ U}+\boldsymbol{u}}(g^{k})-2{\it\alpha}{\it\lambda}\mathscr{M}_{\boldsymbol{ U}+\boldsymbol{u}}(u_{n}^{k}),\end{eqnarray}$$

where $\boldsymbol{u}_{l}^{k+1}$ is the linearized approximation of $\boldsymbol{u}^{k+1}$ . Here, $\mathscr{M}_{\boldsymbol{U}+\boldsymbol{u}}$ is a linearized approximation to $\mathscr{N}$ and $\mathscr{M}_{\boldsymbol{U}+\boldsymbol{u}}(g)$ can be obtained by taking $gf(t,{\it\omega})\boldsymbol{n}$ as the boundary condition, $\boldsymbol{U}+\boldsymbol{u}^{k}$ as the base flow and integrating the linearized NS equation

(2.22) $$\begin{eqnarray}\partial _{t}\boldsymbol{u}-L(\boldsymbol{u})=0\end{eqnarray}$$

from $t=0$ to $t={\it\tau}$ , with the same boundary condition and initial condition as the decomposed NS equation; $\mathscr{M}_{\boldsymbol{U}+\boldsymbol{u}}(u_{n}^{k})$ can be obtained similarly.

Substituting (2.21) into (2.7), we see that the unsteadiness $A$ can be expressed as a polynomial function of ${\it\alpha}$ , and therefore the optimal value of ${\it\alpha}$ (in the linear sense) can be calculated. This optimal value of ${\it\alpha}$ should be substituted into (2.20) to calculate $\boldsymbol{u}^{k+1}$ . When ${\it\alpha}$ is small enough the linearization assumption is satisfied and $A(\boldsymbol{u}^{k+1})<A(\boldsymbol{u}^{k})$ holds. Otherwise, if $A(\boldsymbol{u}^{k+1})>A(\boldsymbol{u}^{k})$ , this linearly optimal value of ${\it\alpha}$ should be reduced.

Since the governing equation (2.5) is nonlinear, the optimal solution of $u_{n}$ calculated above is denoted as the nonlinear optimal control. We see that even when the linearly optimal value of ${\it\alpha}$ is a good estimate of the step length, in each optimization iteration, the decomposed NS equation (2.5) is called once, the linearized governing equation (2.22) is integrated twice and the adjoint equation is called once. If the linear optimal value of ${\it\alpha}$ does not result in a reduced value of the unsteadiness $A$ , then this step length should be reduced and an extra call of the decomposed NS equation is required.

The optimization procedure to calculate the nonlinear optimal control is given in appendix B.

2.4. Sensitivity

We note that when the control cost $E$ is small enough, the adjoint equations can be simplified as

(2.23) $$\begin{eqnarray}\partial _{t}\boldsymbol{u}^{\ast }+L^{\ast }(\boldsymbol{u}^{\ast })+2W^{2}(\boldsymbol{U}-\boldsymbol{U}_{a})=0,\end{eqnarray}$$

indicating that the adjoint variables are functionals of the base flow but independent of the control perturbation, and therefore $g$ is also independent of $u_{n}$ .

Imagining that we start from $u_{n}=0$ , then the solution after one iteration is parallel with $g$ . When updating the solution in the second step following the direction of the gradient, which is also parallel with $g$ , the control does not change, indicating that the optimal control is parallel with $g$ . Therefore, the optimal control can be calculated as

(2.24) $$\begin{eqnarray}u_{n}=-(E/[g,g])^{1/2}g,\end{eqnarray}$$

where the simplified adjoint equation (2.23) is used to calculate the adjoint variables and then the adjoint velocity and pressure are submitted to (2.17) to calculate $g$ . The scalar to scale $g$ in (2.24), i.e.  $-(E/[g,g])^{1/2}$ , is to impose the constraint on the control cost, and the negative sign transforms the control that maximizes unsteadiness to one that instead minimizes unsteadiness. Clearly this calculation only involves a backward integration without forward–backward iterations and therefore the calculation procedure is much simpler than the nonlinear optimal control. This solution is denoted as sensitivity in the following, since it can be interpreted as the sensitivity of the unsteadiness with respect to disturbances on the control boundary. Such an adjoint method has been used to calculate the time sequence of a control signal, which has a fixed spatial distribution and is optimal to suppress linear perturbation growth in a boundary layer flow (Semeraro et al. Reference Semeraro, Pralits, Rowley and Henningson2013).