Michel Bergmann

Laurent Cordier

Jean-Pierre Brancher

Optimal rotary control of the cylinder wake using POD Reduced Order Model

Keywords: numbers: 47.62.+q, 02.60.Pn, 47.27.Vf active flow control, optimal control, Proper Orthogonal Decomposition, Reduced Order Model, Wake flow. 1

In this paper we investigate the optimal control approach for the active control and drag optimization of incompressible viscous flow past circular cylinders. The control function is the time angular velocity of the rotating cylinder. The wake flow is solved in the laminar regime (Re = 200) with a finite element method. Due to the CPU and memory costs related to the optimal control theory, a Proper Orthogonal Decomposition (POD) Reduced Order Model (ROM) is used as the state equation. The key enablers to an accurate and robust POD ROM are the introduction of a time dependent eddy-viscosity estimated for each POD mode as the solution of an auxiliary optimization problem and the use of a snapshot ensemble for POD based on chirp-forced transients.

Since the POD basis represents only velocities, we minimize a drag-related cost functional characteristic of the wake unsteadiness. The optimization problem is solved using Lagrange multipliers to enforce the constraints. 25% of relative drag reduction is found when the Navier-Stokes equations are controlled using an harmonic control function deduced from the optimal solution determined with the POD ROM. Earlier numerical studies concerning mean drag reduction are confirmed: it is shown in particular that without a sufficient penalization of the control input, our approach is energetically inefficient. The main result is that a cost reduction factor of one hundred and 760 is obtained for the CPU time and the memory respectively. Finally, limits of the performance of our approach are discussed.

I. INTRODUCTION

The general aim of active flow control is to alter a natural flow state into another state with more desired properties. Naturally, in order to achieve energetic efficiency of the process, the power needed to control the flow has to stay as low as possible compared to the power saved by the action of control. Flow control has a long history since Prandtl's early experiments for delaying boundary layer separation 1 . However, the recent invention of Micro-Electro-Mechanical Systems 2 , the maturity of control theory and the possible saving of energy 52 that can be offered by an efficient flow control have generated a renewal of interest in active control of fluid dynamical systems 3-5 .

A. Use of approximation models in optimization

Even with recent progress of Computational Fluid Dynamics capability, controlling a fluid flow through a computational approach still remains a formidable endeavor. Due to the non linear character of the Navier-Stokes equations, "real-time" simulation of a fluid system corresponding to a three-dimensional turbulent configuration is unreachable. This situation is even worse in an optimization setting where the large-scale systems (thousands or even millions of degrees of freedom occur frequently in engineering computations) obtained by spatial discretization of the governing equations need to be solved repeatedly. Furthermore, to develop feedback control of complex systems, real-time state solves are necessary 6 . Clearly, reducing the costs of the nonlinear state solutions by some kind of surrogate models 7 of the systems to be controlled is required in flow control. As a result, there have been many studies devoted to the development of Reduced-Order Models (ROM) that serve as lowdimensional approximation models to the large-scale discretized state equations. Some of the reduced-order modelling techniques are restricted to linear systems, other approximation methods can be applied for reducing both linear and nonlinear systems (see Ref. 8 for a review). The model reduction method discussed in this paper fall in the category of reduced basis approaches. This approach consists of seeking an approximation to the state in the form of a linear combination of the so-called reduced basis functions. The coefficients appearing in the linear combination is generally determined through a Galerkin projection of the state equations into the reduced basis space. Ideally, we would like to obtain accurate do so, the reduced basis cannot be determined a priori as it is the case with the finite element method for example. Conversely, the reduced basis is determined a posteriori using experimental or numerical data previously obtained for a given flow configuration, generally for an uncontrolled flow (see Ref. 9 for the consequences in terms of flow control).

For the reduced bases, a number of choices exist (see Ref. 10 for a presentation): Lagrange basis, Hermite basis, Taylor basis, Proper Orthogonal Decomposition (POD) basis, Centroidal Voronoi Tessellations (CVT) basis 11 , etc. Today, the most popular reduced-order modelling approach for complex systems in fluid mechanics is based on POD. Therefore, we will restrict our study to this case and consider that the unsteady non-linear dynamics of the flow is modelled via a reduced order model based on POD (POD ROM). Naturally, the generation of this POD ROM involved an additional numerical cost. First, we need to collect the snapshot ensemble used to determine the POD basis, then we have to solve the corresponding eigenvalue problem and finally to estimate the POD ROM coefficients. Mainly two approaches exist to determine these coefficients. The first, the traditional approach, called POD Galerkin, consists in numerically calculating the coefficients of the system starting from their analytical expression obtained by Galerkin projection. The second, introduced more recently, consists in identifying whole or part of the coefficients of the POD ROM as solutions of minimization problems, one then speaks of calibration procedure or calibrated models (see Ref. 12 for an example). Our hope is that this perhaps large off-line cost can be amortized over several optimization calculations. If the reduced-order modelling is really effective for the controlled flow i.e. if the CPU time necessary to solve the reduced-order model is negligible (of the order of the percent) compared to the corresponding time necessary to solve the high-dimensional model, then optimizing a POD ROM should lead to an important reduction of the computational costs. Moreover, in order to reduce as much as possible this overhead cost, we are interested to use as few runs of the high-dimensional simulation code as possible (ideally one run) in order to generate the snapshots.

The POD (and other similar techniques of ROM) can be viewed as a method of information compression. Essentially, the POD algorithm try to remove "redundant" information (if any) from the data base. As a consequence, the ability of POD modes to approximate any state of a complex system is totally dependent of the information originally contained in the snapshot set used to generate the POD functions. Thus, a POD basis cannot con-tain more information than that contained in the snapshot set. The generation of "good" snapshot set is then crucial to the success of use of POD ROM approach in a bifurcation analysis or more generally in an optimization setting. Since the POD basis is intrinsic to a particular flow, we need to give special attention to adapt the POD ROM (and the POD basis naturally) to changes in physics when the flow is altered by control (control parameters or control function). This central question is illustrated in Fig. 1. Figure 1(a) represents the general configuration of a bifurcation analysis or an optimization problem defined in the control parameter space. In the bifurcation setting, if we knew exactly the bifurcated states then an intelligent sampling of parameter space would consist to take snapshots corresponding to these new operating conditions. This is exactly the method which was used by Ma and Karniadakis 13 to derive a POD Galerkin model able to capture the three-dimensional bifurcation of the wake flow. An equivalent approach consists in incorporating in the POD basis additional vectors called 'shift modes' which point in preferred directions of the phase space. For example, it was shown in Noack et al. 14 that the inclusion of shift modes significantly improves the resolution of the transient dynamics from the onset of vortex shedding to the periodic von Kármán vortex street. In the same paper, it was demonstrated that the inclusion of stability eigenmodes further enhances the accuracy of fluctuation dynamics.

Clearly, some a priori knowledge about the states likely to be represented by the approximation model improves the robustness of the POD ROM. Without any additional information of the different states to be modelled by the POD ROM, the snapshot ensemble for POD should be uniformly distributed in the control parameter space. Then, a lot of runs of the high-dimensional code would be necessary to generate the snapshots and specific methods like the Sequential Proper Orthogonal Decomposition (SPOD) introduced in Ref. 15 would be necessary to develop an accurate model of a controlled transitional flow. In the optimization setting, the problem is almost the same. If we knew at least approximately the path to the optimizer then an intelligent sampling of the parameter space would consist to take snapshots along the path (see Fig. 1(b)). In this interpolatory setting, i.e., if the optimal solution and the path to the optimal solution can be well approximated in the space spanned by the POD basis functions, then it is clear that a POD ROM approach should work. Of course neither the optimized parameters nor the optimal path are known in advance. Then, if we don't have any additional information on the optimal path, there is a high probability that the sampling in parameter space would be unsuitable to approximate correctly the different controlled states encountered by the flow along the optimal path (see Figs. 1(c) and 1(d)). As a consequence, two strategies are conceivable for use of POD ROM in an optimization setting. For lack of systematic, rational, justifiable and effective methodologies for generating good snapshot set (recent work 11 seems to demonstrate that Centroidal Voronoi Tessellations could be one method of intelligent sampling in parameter space), a first approach consists of generating generalized POD functions by forcing the flow with an ad-hoc time-dependent excitation that is rich in transients. The corresponding POD ROM optimization method is schematically described in Fig. 2 where γ(t) denotes the control function defined for the cylinder wake in Sec. II A. The second approach consists of an adaptive method in which new snapshots are regularly determined during the optimization process when the effectiveness of the existing POD ROM to represent accurately the controlled flow is considered to be insufficient. The main drawback of this second approach is that for adaptively updating a reduced basis during an optimization process, new solves of the high-dimensional approximations of the Navier-Stokes equations need to be done. Since these simulations are costly, this approach is not appropriate for real-time control flow. The adaptive method is illustrated in Fig. 3.

B. A prototype of separated flow: cylinder wake

Due to its simple geometry and its representative behavior of separated flows 16 , the viscous flow past a circular cylinder has been extensively used in the past decade as a test-bed to develop methodologies that can be used later to control more complex flows. Different experimental or numerical approaches have been successfully employed for the control of a wake flow but recently optimal control theory attracted increased attention in flow control setting 5,17,18 . For example, He et al. 19 , Homescu et al. 20 , Protas and Styczek 21 used the optimal control theory with the two-dimensional Navier-Stokes equations as the state equation to control by rotary oscillation the unsteady wake of the cylinder (see table I for the characteristics of these approaches). An attractive element of the optimal control approach is the introduction of a cost functional which provides a quantitative measure of the desired objective. However the numerical costs (CPU and memory) associated with the adjoint equation-based methods used to solve these optimization problems are so important that the three-dimensional Navier-Stokes equations are rarely studied 53 . For cutting down the numerical costs different approaches are possible (see Gunzburger 22 for a review). One promising approach is to first develop POD ROM to approximate the fluid flow and then to optimize exactly the reduced models as it was already discussed in Sec. I A. A general discussion of the use of approximation models in optimization can be found in Ref. 23. In this study, we want to develop a low-cost optimal control approach for drag minimization of the cylinder wake with rotary motion for control function (see Fig. 4). Then, to reduce as much as possible the computational costs associated to the present study, the flow is considered two-dimensional and in the laminar regime. However, the methodology presented here that consists of combining the optimal control approach and a POD ROM should easily be expanded to three-dimensional and turbulent flows.

This investigation of drag reduction by unsteady rotary oscillation of the cylinder was motivated in part by the experiments of Tokumaru and Dimotakis 24 where 80% of relative 54 mean drag reduction was empirically found (Re = 15, 000). Recently, Protas and Wesfreid 25 argued (see Sec. II B for numerical evidence and more explanations) that in the supercritical regime of the wake flow, the effectiveness of the control in terms of drag reduction increases with the Reynolds number. Therefore, since the wake flow remains two-dimensional up to a value of the Reynolds number approximately equal to 190 where a spanwise supercritical Hopf bifurcation occurs and where the three-dimensional effects appear 26,27 , the "optimal" value of the Reynolds number for our two-dimensional study is slightly lower than 200.

However for facilitating the comparisons with the results of the literature, a Reynolds number of 200 is considered. According to the observations of He et al. 19 , the control minimizing the drag generates vortices that are less energetic than those produced by the stationary cylinder.

An energetic criterion seems to be well adapted to the investigation of drag reduction.

Therefore, due to the energetic optimality of convergence of the POD basis 28,29 , the choice of POD to develop a reduced order model of the controlled unsteady flow seems to be well adapted. A similar approach was already considered in Graham et al. 30, 31 to control the wake flow at a supercritical Reynolds number of 100.

Finally, we need to choose between the two opposite strategies discussed at the end of Sec. I A. If we want to develop active flow control method that can be used for real-time, on-line feedback control, our interest is to include in the snapshot set all the information needed during the optimization process or at least as much information as we can, and then to generate the reduced order basis. Following this approach the POD functions are realized. This is the specific method that will be used in this paper. The adaptive method has been already tested for the same flow configuration 32 and the results will be published in a subsequent paper.

The main objective of this paper is to emphasize the computational savings of POD ROM based optimal control with respect to more "classical" Navier-Stokes based optimal control as it was already developed for the wake flow in the literature 19-21 . So in this study our main concern is not to determine the control law with the maximum energetic efficiency as it can be characterized for example by the Power Saving Ratio (PSR) (see Protas and Styczek 21 for a definition or hereafter in Sec. VI D) but rather to demonstrate that coupling an optimal control approach and a POD ROM can be successful at least to determine the "optimal" solution corresponding to the open-loop control approach. As far as we know (see table I), the work presented in Ref. 21 is the only one which considers for cost functional the sum of the drag power and the control power thus making it possible to determine an optimal solution that is by construction energetically efficient. In the other works, the cost of the control is not considered or at best as a regularization parameter. This discussion will be developed in Sec. VI D where we compare the energetic efficiency of the different approaches. This article is organized as follows: Sec. II introduces the flow configuration and describes the numerical method used to simulate the flow. In the next two sections, the Proper Orthogonal Decomposition is first introduced (Sec. III), then we outline the control function method used to develop a POD ROM of the controlled flow (Sec. IV). The optimal control problem is stated in Sec. V which includes the definition of the cost functional and a description of the Lagrange multiplier method used to solve the constrained optimization problem. Finally, before to present the results of the POD ROM based control in Sec. VI B and the drag reduction obtained with the Navier-Stokes equations when the optimal control function determined with the POD ROM is used (Sec. VI C), we describe in Sec. VI A how to determine generalized POD basis functions.

II. PROBLEM FORMULATION AND SIMULATION METHOD

A. Flow configuration, governing equations and numerical method

Let Ω be a two-dimensional bounded region filled with a Newtonian incompressible viscous fluid of kinematic viscosity ν. We denote by Γ the boundaries of Ω and we note U ∞ the uniform velocity of the incoming flow (see Fig. 4). Wake flows dynamics are characterized 16 by the Reynolds number Re and by the natural Strouhal number St n at which vortices are shed in the wake of the cylinder. Traditionally, the Reynolds number is defined as Re = U ∞ D/ν where D is the cylinder diameter (R is the corresponding radius) and the natural Strouhal number as St n = f D/U ∞ where f is the fundamental frequency. The rotary control is characterized by the instantaneous rotation rate θ(t). Equivalently, we can specify the tangential boundary velocity V T (t) = R θ(t) or the non dimensional velocity γ(t) defined as the ratio of the tangential velocity V T (t) to the upstream velocity U ∞ where, here, for numerical convenience, U ∞ is supposed to be equal to unity. Hereafter, the control function γ(t) is sought using the optimal control theory in order to minimize the mean drag coefficient of the wake flow.

The flow of viscous incompressible fluid is governed by the Navier-Stokes system that can be expressed (u ⊗ v denotes the dyadic product of the two vectors u and v, i.e., (u ⊗ v) ij = u i v j ) in non dimensional 55 form as

∇ • u = 0 in Ω × (0, T ), ∂u ∂t + ∇ • (u ⊗ u) = -∇p + 1 Re ∆u in Ω × (0, T ), (1)
with an initial condition defined as

u| t=0 = u 0 with ∇ • u 0 = 0, (2)
where u = {u i (x, y, t)} 2 i=1 is the two-components fluid velocity vector, p = p (x, y, t) the pressure field and t the time. In the above equations (0, T ) is the time interval during which the flow is considered.

The problem specification is now completed by the boundary conditions. At the left boundary, an inflow boundary condition is applied:

(u 1 , u 2 ) = (1, 0) on Γ i × (0, T ).
(3)

At the transverse boundaries, zero shear stress conditions are imposed:

∂u 1 ∂y = 0, u 2 = 0 on Γ b × (0, T ) and Γ t × (0, T ). (4)
At the outflow boundary, a non-reflecting boundary condition is considered. The velocity field is deduced on Γ o as the solution of a wave-like equation 32 :

∂u ∂t + u 1 ∂u ∂x - 1 Re ∂ 2 u ∂y 2 = 0 on Γ o × (0, T ). (5)
When this artificial boundary condition is used, no spurious reflections from the downstream boundary are observed thus making it possible to reduce the computational domain in the cylinder wake.

Finally, on the cylinder surface the velocity is equal to the tangential boundary velocity.

Since

U ∞ = 1, this condition becomes u(x, t) = γ(t)e θ (x) on Γ c × (0, T ), (6)
where e θ is the unit tangent vector on Γ c .

The partial differential equations ( 1) are discretized in time by a three steps projection method and in space using a Galerkin finite element approximation (P 1 , P 1 ). This numerical method is classic and the details that can be found in Cordier and Bergmann The objective of this section is to demonstrate that for the laminar regime considered in this paper, the two-dimensional numerical simulation discussed previously can be viewed to represent correctly the dynamics of the cylinder wake flow. The well-known Von Kármán street, characteristic of the cylinder wake, is clearly visible in Fig. 6 In a viscous flow the total forces acting on a body are contributed by the pressure and skin friction terms. Let K p be the pressure coefficient defined by

K p = 2(p -p ∞ )
where the subscript ∞ denotes quantities evaluated on the inflow boundary Γ i . The aerodynamic coefficients can then be calculated as:

C(t) = - Γc K p n dΓ + 2 Re Γc ∂u ∂n dΓ = C D (t) e x + C L (t) e y (
C D = C base D + C 0 D . (8)
As it can be viewed in Fig. 9, when the Reynolds number increases, the relative contribution of the unsteady part of the flow to drag becomes more significant. At a given Reynolds number, the contribution of the base flow to drag cannot be modified. Then controlling the wake flow by rotary oscillation can only reduced the contribution of the unsteady part.

Since the contribution C 0 D increases with the Reynolds number, the controllability of the flow is more important for higher Reynolds numbers as it was announced in Sec. I B. This discussion points out a fundamental question already arisen in Ref. 25: is it possible to assert that the minimal value of drag that can be obtained under periodic forcing conditions and, at a given value of the Reynolds number, corresponds to the base flow solution? In other words, is it possible to obtain a mean flow correction with negative drag? Certainly, this question deserves more study from the computational and theoretical points of views.

However, we report in Bergmann 32 that when the rotation is applied only on a limited part of the cylinder boundary, a mean drag coefficient lower than the corresponding base flow solution was obtained.

Finally, in Fig. 10 we represent the streamlines of the base flow obtained at Re = 200. In order to determine this unstable flow with (almost) the same numerical code used to obtain the unsteady flows, we impose symmetry conditions in the wake at every time step of the simulation, inhibiting in such a way the growth of any symmetry-breaking perturbations. This kind of method was already used in Ref. 25 for the same purpose.

III. PROPER ORTHOGONAL DECOMPOSITION (POD)

This method was introduced in turbulence by Lumley in 1967 as an unbiased definition of the coherent structures widely known to exist in a turbulent flow. Starting with a set of realizations of the velocity fields u(X) where X = (x, t) ∈ D = Ω × R + , a coherent structure is defined as the deterministic function Φ(X) which is most similar on average to the realizations u(X). Mathematically, the notion of most similar corresponds to the solution of the following constrained maximization problem:

max Φ |(u, Φ)| 2 subject to Φ 2 = 1 (9)
where (., .) denotes a scalar product in the Hilbert space of square-integrable functions L 2 , . is the corresponding norm and the brackets . denote an averaging operation, which may be a time or space average. More details on POD and all the justifications can be found in Cordier and Bergmann 29 .

From variational calculus it can be shown that the problem ( 9) is equivalent to a Fredholm integral eigenvalue problem:

D R ij (X, X ′ )Φ j (X ′ ) dX ′ = λΦ i (X) (10)
where R ij (X, X ′ ) is the two-point space-time correlation tensor (here and in the following, i and j vary from 1 to n c where n c is the number of velocity components). Since R ij is self-adjoint and non-negative definite, it follows from the Hilbert-Schmidt theory that equation ( 10) has a denumerable infinite number of eigenvalues λ n and eigenfunctions Φ i n (n = 1, . . . , +∞). These eigenvalues are all real and positive and form a decreasing and convergent series. Each eigenvalue represents the contribution of the corresponding modes Φ n to the total kinetic energy. The associated eigenvectors Φ n (also called empirical eigenfunctions) form a complete orthogonal set and have been normalized, so that they verify

(Φ n , Φ m ) = δ nm .
Moreover, it can be demonstrated 29 that they are optimal in an energetic sense (for a given number of modes N , the projection on the subspace spanned by the N leading eigenfunctions will contain the greatest possible kinetic energy on average).

Depending on the choice made for the average operator . appearing in ( 9), two equivalent formulations of POD can be found 29 . When the average is estimated in time, the first approach called classical POD or direct method and originally introduced by Lumley is obtained. In this case, the kernel R ij of the Fredholm equation ( 10) is replaced by the two-point spatial correlation tensor r ij (x, x ′ ) and the eigenfunctions Φ(X) by φ(x). In the second case suggested by Sirovich 35 and called snapshot POD, the average operator is evaluated as a space average over the domain in interest. The Fredholm equation to be

T C(t, t ′ )a n (t ′ ) dt ′ = λ n a n (t) (11)
where C(t, t ′ ) is the temporal correlation tensor constructed as

C(t, t ′ ) = 1 T (u(x, t), u(x, t ′ )) (12)
where the outer parentheses (., .) represent the inner product defined as

(u, v) = Ω (u, v) 2 dx = Ω nc i=1 u i v i dx.
In Eq. ( 11), a n are the time-dependent POD eigenfunction of order n. These modes form an orthogonal set, satisfying the condition:

1 T T a n (t)a m (t) dt = λ n δ nm . ( 13
)
The spatial basis functions φ i n can then be calculated from the velocities u i and the coefficients a n with:

φ i n (x) = 1 T λ n T u i (x, t)a n (t) dt. (14)
Since the POD eigenfunctions can be represented as linear combinations of the velocity fields, they inherit all the properties of the original data that are linear or homogeneous.

Hence the eigenfunctions are divergence free for an incompressible fluid (∇ • φ n = 0) and verify automatically the homogeneous boundary conditions of the numerical simulation used to determine the flow realizations.

For reasons of statistical convergence of the average operator, the snapshot POD is more appropriate when data issued from numerical simulations are used. Hence, this method was adopted in this work.

Finally the set of POD modes {φ n } +∞ n=1 is complete in the sense that any velocity field u(x, t) can be expanded in the eigenfunctions as

u i (x, t) = N P OD n=1 a n (t)φ i n (x) (15)
where N P OD is equal to the number of flow realizations used to solve the POD problem (11).

IV. POD ROM OF THE CONTROLLED CYLINDER WAKE

When the rotary control is applied, the boundary conditions on the cylinder become inhomogeneous and time-dependent. As a consequence, the POD basis functions used in the Galerkin projection have not homogeneous boundary conditions and extra terms appear in the POD ROM (see the pressure term in Eq. 21). To address this situation, the control function method introduced in Graham et al. 30 is used.

A. The control function method

Here, we decide to justify the introduction of a suitable "control function" u c to remove the inhomogeneous boundary conditions on Γ in the general context of boundary control problem for fluid flows 33 . So, we consider the Navier-Stokes system (1) completed with the Dirichlet boundary conditions:

u(x, t) = h(x, t; U (t)) on Γ × (0, T ), ( 16
)
where U is the control input.

Finally, we assume that the boundary of the domain, Γ, can be split into two parts such that Γ c denotes that part of the boundary where the control is applied and Γ \ Γ c is the part of the boundary that is not controlled. More precisely, in the case of the controlled cylinder wake flow that is considered, the boundary conditions can be written as (see Eq. 6):

h(x, t; U (t)) =      γ(t)e θ (x) on Γ c × (0, T ), c(x) on Γ \ Γ c × (0, T ). (17)
Since the boundary conditions h depend on time, the velocity expansion is now defined as

u(x, t) = u m (x) + γ(t) u c (x) + u(x, t) (18) with u(x, t) = N P OD k=1 a k (t)φ k (x) (19)
where u m (x) is the mean velocity field obtained as an ensemble average of the flow realizations contained in the snapshot set and where u c (x) is a reference flow field, called control function, that describes how the control action γ(t)e θ (x) influences the flow. This control function satisfies the following boundary conditions:

γ(t)u c (x) =      γ(t)e θ (x), on Γ c × (0, T ), 0 on Γ \ Γ c × (0, T ). (20)
A convenient way to generate it is to take the solution of the governing equations (1) for the steady cylinder rotation corresponding to γ = 1 and homogeneous boundary conditions for the uncontrolled part of Γ (see Fig. 11 for an illustration).

When the snapshot POD approach is considered for deriving a POD based reduced order model, the input ensemble used to determine the POD modes consists of N t flow realizations called time snapshots u(x, t i ), x ∈ Ω, taken at time instants

t i ∈ (0, T ) , i = 1, • • • , N t .
In case of time-dependent boundary conditions, the procedure for computing the POD basis can be formulated as follows 33 . A mean velocity u m (x) is first computed as the ensemble average of the modified input data defined as

U ′ = {u(x, t 1 ) -γ(t 1 )u c (x), • • • , u(x, t Nt ) - γ(t Nt )u c (x)}.
Afterward, the POD basis functions φ k are estimated with the input collection boundary, the contribution of the pressure term in Γ o is not exactly zero for the cylinder wake. However, for this flow configuration, the Galerkin projection of the pressure term is found to be negligible (see Noack et al. 14 and the discussion in Appendix A 1).

U ′′ = {u(x, t 1 ) -γ(t 1 )u c (x) -u m (x), • • • , u(x, t Nt ) -γ(t Nt )u c (x) -u m (x)}. Since (u(x, t i ) -γ(t i )u c (x))| Γc =

B. Derivation of the POD ROM

The weak form of the Navier-Stokes equations is restricted to the POD subspace S P OD N gal spanned by the first N gal spatial eigenfunctions φ i . The energetic optimality of the POD basis functions suggests that only a very small number of POD modes may be necessary to describe efficiently any flow realizations of the input data. The dimension N gal ≪ N P OD of the subspace S P OD N gal is the smallest integer M such that the Relative Information Content (RIC) defined as the ratio M i=1 λ i / N P OD i=1 λ i is greater than δ% where δ is a predefined percentage of energy (here δ = 99 and M = N gal = 4 for the uncontrolled flow, see Fig. 18).

The Galerkin projection yields 33 :

φ i , ∂u ∂t + (u • ∇)u = (p, ∇ • φ i ) -[p φ i ] - 1 Re ∇φ i , (∇u) T + 1 Re (∇u) T φ i (21) with [u] = Γ u • n dx and A, B = Ω A : B dx = nc i, j=1 Ω A ij B ji dx.
Inserting the expansion (18) into the Galerkin projection ( 21) of the Navier-Stokes equations, we obtain after some algebraic manipulations the reduced order control model:

d a i (t) d t =A i + N gal j=1 B ij a j (t) + N gal j=1 N gal k=1 C ijk a j (t)a k (t) + D i d γ d t +   E i + N gal j=1 F ij a j (t)   γ + G i γ 2 i = 1, • • • , N gal . (22)
The coefficients

A i , B ij , C ijk , D i , E i , F ij and G i depend explicitly on φ, u m and u c .
Their expressions are given in Appendix A.

The system of equations ( 22) is then integrated in time with a fourth order Runge-Kutta scheme from a given set of initial conditions

a i (0) = (u(x, 0), φ i (x)), i = 1, • • • , N gal (23)
yielding a set of predicted time histories for the mode amplitudes a i (t) which can be compared with the POD temporal eigenfunctions.

Due to the truncation involved in the POD-Galerkin approach, the higher POD modes corresponding to the dissipative scales of the flow are not explicitly taken into account in the POD ROM. As a consequence, when the equations ( 22 and to improve the accuracy of POD flow models 38 . Here, the low-dimensional Galerkin model ( 22) is stabilized by a time dependent eddy-viscosity estimated for each POD mode as the solution of an auxiliary optimization problem described in Bergmann 32 .

As shown in Fig. 12 for an uncontrolled flow (γ = 0), when the POD ROM is stabilized numerically, excellent qualitative and quantitative agreements are found between the integrated time histories of the POD modes kept in the truncation (predicted modes) and the results obtained by the numerical simulation (projected modes). For a controlled flow, an accurate description of the dynamical behavior is still possible (see Fig. 13). However, if the objective is to use the POD ROM in an optimization process, special care is needed to collect the set of snapshots used to determine the reduced basis (see the discussion in Sec. I A). This point is particularly described in Sec. VI.

V. OPTIMAL CONTROL APPROACH

In this section we discuss how the optimal control approach can be used to determine the rotation rate γ(t). The aim is to minimize a cost functional J , which incorporates the control goal and some measure of the control effort, over a certain period of time T o corresponding to few periods of the Von Kármán street (T o is generally referenced as the optimization horizon). Here we envisage employing the POD ROM of Sec. IV for model-based control of the vortex shedding flow. Therefore, since only the flow velocities are directly represented by the POD basis functions and since the pressure drag represents about 80% of the total drag, our objective is to minimize a drag-related cost function. A natural control aim is the reduction of the wake unsteadiness i.e. the energy contained in the wake as defined in Ref. 31. Mathematically, this goal is expressed as the following functional

J ( u, γ(t)) = To 0 Ω J( u(x, t), γ(t)) dxdt = ℓ 1 2 To 0 Ω u(x, t) 2 2 dxdt + ℓ 2 2 To 0 γ 2 (t) dt
where the first term represents the control goal and the second a penalization term. In this formulation ℓ 1 and ℓ 2 are two positive regularization parameters that can be empirically chosen to limit the size of the control. This formulation is equivalent to the introduction in front of the cost control term of a tuning parameter ℓ defined as the ratio of ℓ 2 to ℓ 1 as can be found for example in Ref. 39. If ℓ is chosen to be large, the control is "expensive", if it is small, the control is "cheap". In other words, the parameter ℓ weighs the control cost with respect to the state energy.

Hereafter, let a be the vector containing the first N gal time-dependent expansion coeffi-

cients a i (a = a 1 , a 2 , • • • , a N gal ).
Introducing the POD expansion ( 19), the functional J becomes:

J (a, γ(t)) = To 0 J(a, γ(t)) dt = ℓ 1 2 To 0 N gal i=1 a i 2 (t) dt + ℓ 2 2 To 0 γ 2 (t) dt. (24)
The flow control problem is then expressed as:

         min γ(t) J (a, γ(t)) subject to N (a, γ(t)) = 0 (25)
where the constraints N (a, γ(t)) = 0 correspond to the POD ROM ( 22).

The constrained optimization problem (25) is solved using the Lagrange multiplier method as described in Gunzburger 40 . The constraints are enforced by introducing the Lagrange multipliers or adjoint variables ξ and the Lagrangian functional

L(a, γ, ξ) = J (a, γ(t)) -(ξ, N (a, γ)) 2 t (26) = J (a, γ(t)) - N gal i=1 To 0 ξ i (t) N i (a, γ) dt.
The solutions (states a, co-states ξ and control γ) of this new unconstrained optimization problem are such that L(a, γ, ξ) is rendered stationary:

δL = N gal i=1 ∂L ∂a i δa i + ∂L ∂γ δγ + N gal i=1 ∂L ∂ξ i δξ i ≡ 0
where δa, δγ and δξ are arbitrary variations.

Considering 56 that each argument of L is independent of the others, the optimality system is determined by setting the first variation of L with respect to ξ, a and to γ to be equal to zero.

Setting the first variation of L with respect to the Lagrange multiplier ξ equal to zero, we recover the state equation N (a, γ(t)) = 0.

We now set to zero the first variation of L with respect to the state variable a. After integration by parts, the adjoint equations

d ξ i (t) dt = -ℓ 1 a i (t) - N gal j=1   B ji + γ(t) F ji + N gal k=1 (C jik + C jki ) a k (t)   ξ j (t) i = 1, • • • , N gal (27)
with the terminal conditions:

ξ i (T o ) = 0, i = 1, • • • , N gal (28)
are derived.

Finally, setting the first variation of L with respect to the control γ equal to zero yields the optimality conditions

δγ(t) = - N gal i=1 D i dξ i dt + ℓ 2 γ + N gal i=1   E i + N gal j=1 F ij a j + 2G i γ(t)   ξ i . (29)
Since the optimal control is the control that minimizes J , i.e. that nullifies the gradient of the cost functional with respect to the control variable ∇ γ J , it can be shown 40 that imposing ∇ γ J = 0 is equivalent to solving the optimality conditions (29), therefore

∇ γ J (t) = δγ(t).
The first-order necessary conditions yield a system of coupled ordinary differential equations (state equation 22, adjoint equations 27 and optimality condition 29) called optimality system. Due to large storage and CPU costs that system is rarely solved without iteration. Instead a simple iterative process can be effected as follows (in this algorithm and in Fig. 14, the superscripts (n) denote the iteration number):

Start with an initial guess γ (0) (t) for the control function γ(t). For n = 0, 1, 2, • • • , 1. Solve the POD ROM ( 22) forward in time to obtain the corresponding mode amplitudes a (n) (t).

2. Use a (n) (t) computed in step 1 to solve the adjoint equations ( 27) backward in time for the adjoint variables ξ (n) (t).

3. Use the state variables a (n) (t) computed in step 1 and the adjoint variables ξ (n) (t) computed in step 2 to estimate the optimality condition ( 29) and determine the functional

gradient ∇ γ J (n) (t) on the interval [0; T o ].
4. Use this estimation of gradient to update the control γ (n+1

) (t) = γ (n) (t) + α (n) d (n) (t).
Here d (n) is a direction of descent estimated with one's favorite optimization method using the gradient of the functional ∇ γ J (n) (t) and α (n) is the length step in that direction.

5. If some stopping criterion is satisfied, stop; otherwise, return to step 1.

VI. POD ROM BASED CONTROL

In Sec. IV, it was demonstrated that at least at the design conditions the stabilized POD ROM (22) represents accurately the dynamics of the controlled flow. Therefore in this section, we use the results of the optimal control theory presented in Sec. V to determine a model based control function γ(t).

A. Generalized POD basis functions

In this paper we make the choice to not refresh the POD ROM ( 22) during the optimization process (see Sec. I B). Clearly, it corresponds to the assumption that the path to the optimal solution and the optimal solution can be well approximated by the original POD basis, hence the importance of generating "good" snapshot sets. Here, we follow the method introduced in Graham et al. 30 and derived generalized POD basis functions that correspond to an ad-hoc forcing term rich in transients. More exactly, we impose as rotation rate for the cylinder a slowly varying amplitude and frequency sinusoid or 'chirp' function. This temporal excitation γ e shown in Fig. 15 is mathematically represented by the function:

γ e (t) = A 1 sin(2πSt 1 t) × sin(2πSt 2 t -A 2 sin(2πSt 3 t))
where A 1 = 4, A 2 = 18, St 1 = 1/120, St 2 = 1/3 and St 3 = 1/60.

Although there is some arbitrariness involved in the choice of the parameters in the definition of γ e , the use of a snapshot ensemble based on chirp-forced transient represents an attractive compromise between the number of high-dimensional runs necessary to generate the snapshots (one) and the robustness of the POD basis. However, we can hardly imagine that the generalized POD basis functions derived with these values of the parameters can accurately reproduce the controlled dynamics corresponding to an important penalization of the control input (parameter ℓ 2 in Eq. 24) i.e. to low amplitudes of the forcing (for comparison in Ref. 21, where the sum of the work needed to resist the drag force and the work needed to control the flow was considered as cost functional, the optimal amplitude A is always found lower than 0.5). This behavior will be illustrated in Sec. VI B.

In Fig. 16 The control function γ is determined as a converged solution of the iterative process introduced at the end of Sec. V. As an initial guess for the control function the excitation γ e is selected (see Fig. 14). In step 4 of the iterative method the Fletcher-Reeves version of the Conjugate Gradient Method (see Ref. 41 for example) is used. In every iteration the control function is updated according to

γ (n+1) (t) = γ (n) (t) + α (n) d (n) (t), (30)
where d (n) represents the conjugate direction given by

d (n) (t) = -∇ γ J (n) (t) + β n d (n-1) (t), d (0) (t) = -∇ γ J (0) (t) (31)
with β n a coefficient given by

β n = ∇ γ J (n) , ∇ γ J (n) (∇ γ J (n-1) , ∇ γ J (n-1) ) . ( 32
)
The linear search parameter α (n) is computed at each iteration by the backtracking Armijo method 42 , an algorithm that assures that the corresponding step is not too small and verifies the Goldstein condition. The iterative method is stopped when two following values of the functional J are sufficiently close i.e. when |∆J (a, γ)| = |J (n+1) (a, γ) -J (n) (a, γ)| < 10 -5 .

Once this criterion of convergence is reached, the relative reduction of the cost functional J characterizing the wake unsteadiness is equal 57 to 69% for ℓ 2 = 0. Figure 20 represents the variation of the corresponding cost functional J with respect to the iteration number n during the optimization process. The decrease of the cost functional is very rapid: though a number of 600 iterations is necessary to obtain the above mentioned value of relative reduction, 60% of decrease is already obtained for the first 30 iterations. The time evolution of the corresponding wake unsteadiness as defined by Eq. ( 24) is represented in Fig. 21. An important decrease of the wake unsteadiness is observed when the optimal control function γ opt is applied. The same observation can be realized by comparing the time evolution of the first six POD modes predicted by the POD ROM ( 22) when the control function is equal to γ e (Fig. 22) and when the optimal control function γ opt is used (Fig. 23). Clearly, the amplitudes of the temporal POD functions are considerably reduced. Finally the time evolution of the optimized control function γ opt corresponding to ℓ 2 = 0 is shown in Fig. 24.

optimally controlled flow (γ = γ opt ) when ℓ 2 = 0. The relative mean drag reduction was found to be of the order of 25% (from approximately an average value of 1.4 to an average value of 1.04). This reduction of mean drag coefficient is substantial but the optimal control parameters A and St f found by the optimization process based on the POD ROM and the wake unsteadiness for cost functional (A = 2. Here, when ℓ 2 = 0 the optimal control γ opt is able to annihilate about 77% of the drag related to vortex shedding.

As it can be intuitively expected, when ℓ 2 = 1 the relative mean drag reduction is found to be lower than when ℓ 2 = 0. The reduction is of the order of 8%, the mean drag coefficient varies approximatively from a value equal to 1.4 for the uncontrolled wake to a value equal to 1.29 when the flow is controlled with a sinusoidal rotation defined with A = 0.75 and

D. Discussion

The numerical results obtained here with the POD ROM as state equation agree to a large extent to results obtained in other numerical approach where the optimal control theory is applied for the same flow configuration directly to the Navier-Stokes equations (see table I for the characteristics of the different algorithms). However, quantitative comparisons of the control algorithms are difficult because for the comparisons to be fair, it would be necessary that the same actuation method and the same control objectives were used in the various studies. Hereafter, we will see that the studies used for comparison were performed with either a different actuation method, or a different control objective.

Therefore, only qualitative comparisons of the control methodologies are possible.

Protas and Styczek 21 obtain, with a rather small magnitude of the control (γ(t) ≤ 0.2), a drag reduction of about 15% for a Reynolds number equal to 150, presenting a less significant controllability. Contrary to our findings, the time history of their optimal controls exhibit spiky behaviors difficult to represent by an harmonic oscillation. The main reason is the following. Due to the numerical costs related to the use of the Navier-Stokes equations, Protas and Styczek were constrained to use a particular approach referred to as piecewise optimal. In this approach, the optimization interval T o is cut in a sequence of shorter intervals where the optimization is independently performed. Simply, the state reached by the optimized flow at the end of a given interval is taken as the starting point for the optimization in the following interval. Of course, as it was already remarked by these authors, optimal controls found by this approach does not necessarily correspond to the solution of the optimization problem defined on the whole optimization time-span. However, they presented results that indicate that an optimization horizon comparable to the length of the natural vortex shedding period is sufficient to decrease the mean drag. In our case (see Although the relative mean drag reduction obtained by Protas and Styczek 21 is the weakest among all that found by the different authors, we can note that from an energetic point of view it is their approach which is the most effective. Indeed, the drag modifications were obtained at different costs in the different approaches discussed above. Hence, for the comparison of the relative drag reduction to be fair, one must take into account the costs of the control. The energetic efficiency of the control is represented as the ratio of saved power to input power. Following Protas and Styczek 21 , it can be characterized by the Power Saving Ratio defined as:

P SR = P D uncontrolled T -P D controlled T P C T , (33)
where P D T and P C T represent respectively the mean of the instantaneous drag power P D and control power P C estimated over a finite horizon T . In our configuration, the drag power can be easily deduced from the drag coefficient C D since P D = F D U ∞ where F D is the drag force exerted on the cylinder. Moreover, it can be shown (see Protas and Styczek 21 ) that P C = M z θ where M z is the torque applied to the cylinder and θ the angular velocity.

In the present study, the PSR are found equal to 0.25 when ℓ 2 = 0 (of the same order of the value obtained in the open-loop control study of Protas and Wesfreid 25 ) and to 0.86 i.e. still lower than unity when ℓ 2 = 1. So an harmonic rotary control with the control parameters found in our approach is energetically inefficient. In their papers, He et al. 19 and Homescu et al. 20 did not give information on the energetic efficiency of their algorithms. However, we can estimate that with their values of optimal control parameters, their PSR was lower than unity as well. For comparison, Protas and Styczek 21 found a value of PSR equals to 51 at Re = 150. Naturally this energetic efficiency is explained by the choice particularly adapted to the PSR criterion of the total power (drag and control) as cost functional. One way to improve the PSR in our study consists to divide the cylinder boundary into two regions: an upstream part on which an optimal control law is applied, and the remainder of the cylinder which is not controlled (see Bergmann 32 for numerical evidence).

Except for the approach presented in Ref. 21 where the energetic efficiency is favored, the drag reduction found with the Navier-Stokes equations as state equation is only slightly higher than the one found with our approach but the numerical costs (CPU and memory) associated to their control are much more important. Using a POD ROM the costs of the flow solves necessary at each iteration of the optimizer are greatly reduced. In our study, the CPU time necessary to obtain with the POD ROM the flow dynamics over one natural vortex shedding period represents 1% of the time necessary to solve the Navier-Stokes equations with the finite-element approach. In first approximation 58 the same gain is obtained for the adjoint equations and the optimality conditions. The total CPU cost is thus drastically reduced (approximately a factor equal to one hundred). With regard to memory cost, note that we need to store the latest state approximation for all space-time to solve the adjoint equations and all the adjoint variables to estimate the optimality conditions. When the finite element simulation is used to solve the optimal control problem over a time horizon T o , we need to store the state and adjoint variables (two velocity components and the pressure) at every time-step and for each vertex of the mesh. When the POD ROM is used, we only need to store the time evolution of the state variables a and of the adjoint variables ξ for N gal POD modes plus the coefficients appearing in the state equation ( 22) (eventually the POD basis functions can be stored for a future use). The parameters used in Sec. VI B are T o = 30 for the time horizon (approximately six times the natural vortex shedding period), ∆t = 0.01 for the optimization time-step, N v = 13, 000 for the number of vertices and N gal = 40 for the number of POD modes kept in the ROM. After estimation we found that the memory cost of the POD ROM approach is approximately 760 times lesser than for the Navier-Stokes model (approximately 280 if we decide to store the POD eigenfunctions). The reduction of the numerical costs offered by our approach is so important that the study of three-dimensional unsteady complex flows by the optimal control theory becomes possible.

However, as it was suggested by Gunzburger in Ref. 22, the success of our approach depends on the ability of the POD basis to well approximate the optimal solution and the path to the optimal solution. Using a POD ROM to solve an optimization problem in the extrapolary regime is not so clear. Certainly, some updating of the POD basis would be necessary during the optimization process like in the Trust Region POD method introduced by Fahl 45 . We already evaluated this method for the same flow configuration 32 and the results will be published elsewhere.

VII. CONCLUSIONS

The objective of this paper was to illustrate the potential gain that can be offered by the use of the Proper Orthogonal Decomposition for optimal control of fluid flows. Our methodology was presented for the unsteady rotary control of the cylinder wake in the laminar regime (Re = 200). Defining a cost functional representative of the wake unsteadiness, the optimal control problem was solved with a POD ROM of the controlled flow as the state equation. The solution of the optimization process was then used to control numerically the wake flow with the Navier-Stokes equations as flow model. Finally, a significant reduction (25%) of the amplitude of the drag coefficient was found. However, we demonstrate that our approach is energetically inefficient. These numerical results agree to a large extent to results obtained by other researchers 19-21 using the two-dimensional Navier-Stokes equations to solve the optimal control problem. Comparing to those studies, the main advantage of our approach is that the numerical costs (CPU and memory) are negligible (of the order of 1% for the CPU time and even less for the memory cost). The conceptual drawback is that there is no mathematical assurance that the solution of the optimization algorithm working with the approximation models will correspond to the solution of the optimization problem for the original dynamical system. As it was suggested by Alexandrov et al. in Ref. 23, a possible way to be assured that the solution of the optimization problem for the reduced order model is likely to yield at least to a local optimum for the original high-fidelity problem, is to use the general trust region framework 59 . Therefore is the POD ROM approach useful in the flow optimization setting ? A partially answer may be given by quoting

Gunzburger 40 "without an inexpensive method for reducing the costs of flow computations, it is unlikely that the solution of optimization problems involving the three-dimensional, unsteady Navier-Stokes system will become routine".

equal to zero. Hence, the Galerkin projection of the pressure term (A2) writes

(φ i , -∇p) = - Γo p φ i • n dx. ( A3
A i = -(φ i , (u m • ∇)u m ) - 1 Re (∇ ⊗ φ i ) T , ∇ ⊗ u m + 1 Re [(∇ ⊗ u m )φ i ] , B ij = -(φ i , (u m • ∇)φ j ) -(φ i , (φ j • ∇)u m ) - 1 Re (∇ ⊗ φ i ) T , ∇ ⊗ φ j + 1 Re [(∇ ⊗ φ j )φ i ] , C ijk = -(φ i , (φ j • ∇)φ k ) , D i = -(φ i , u c ) , E i = -(φ i , (u c • ∇)u m ) -(φ i , (u m • ∇)u c ) - 1 Re (∇ ⊗ φ i ) T , ∇ ⊗ u c + 1 Re [(∇ ⊗ u c )φ i ] , F ij = -(φ i , (φ j • ∇)u c ) -(φ i , (u c • ∇)φ j ) , G i = -(φ i , (u c • ∇)u c ) .
An exception is the seminal work of Bewley et al. 41 where the optimal control theory is used to determine controls that reduce the drag of a turbulent flow in a three-dimensional plane channel simulated at Re τ = 180.

Let C D (t) be the instantaneous drag coefficient defined in Eq. ( 7) and C D T be the mean temporal drag coefficient estimated over a finite horizon T , the relative mean drag reduction is defined as

C D unforced T -C D forced T C D unforced T ,
where the terms 'unforced' and 'forced' are used respectively for non-rotating and rotating cylinder.

In the following, all variables are assumed to be non-dimensionalized with respect to the cylinder diameter D and the oncoming flow U ∞ .

This was not true for the original optimization problem (25) involving J since the arguments a and γ were constrained to satisfy N (a, γ(t)) = 0.

The value of J is 11.85 at the beginning of the iterative process and 3.70 when convergence is obtained. To be used as a reference value, J is approximately equal to 10 for an uncontrolled flow.

58 Contrary to the Navier-Stokes equations, the adjoint and optimality equations are linear in the adjoint variables and hence the numerical costs necessary to solve these equations are usually lower than for the Navier-Stokes equations. However, since the adjoint equations are defined in reverse time, the adjoint fields may grow exponentially when they are calculated, requiring the use of regularization methods as recently discussed in Ref. 48. In this study, no exponential growth was noticed with the adjoint equations ( 27).

59 The trust region mechanism gives a measure of how well the approximation model is predicting improvement in the high-fidelity model and thus suggests criteria for automatically changing or improving the reduced model when poor or incorrect prediction of improvement is obtained.

a i (t), γ (n) (t)
a i (t), γ (n) (t), ξ i (t)

∇ γ J (n) (t)

  33 will not be discussed here. The discrete equations are numerically solved on an unstructured mesh with the Partial Differential Toolbox of Matlab. The accuracy of the numerical code was extensively tested in Bergmann 32 for different time steps, mesh sizes and Reynolds numbers varying from Re = 4 (creeping flow) to Re = 1000. The dependence of the mean drag coefficient and the natural Strouhal number on the Reynolds number were evaluated by comparison of reference results available in the literature. In particular, the well-known over-prediction of the drag coefficient for two-dimensional simulations 16 was observed at Re = 1000. Numerically, it was found that for a time step equal to ∆t = 1.5 10 -2 and a finite element mesh consisting of 25, 000 triangles and 12, 686 vertices (see Fig.5), the present simulations described accurately the dynamics of the uncontrolled and controlled flows (in appendix B are presented typical results of an open-loop control study of the cylinder wake). Therefore, in the following, we describe only the results of the simulations performed at Re = 200. B. Results of the simulation at Re = 200

  (a). Moreover, since non reflective boundary conditions have been considered on Γ o , no spurious reflections from the downstream boundary occur (see Fig. 6(b)).

  7)where C D and C L represent respectively the drag and lift coefficients and n is the outward unit normal at the boundary surface.The time histories of the lift and drag coefficients after the long term behavior has been established are represented in Fig.7. The time mean value of C D is 1.39 and the unsteady amplitudes of C L and C D are 0.0921 and 1.38 respectively (see tableII). The periodic regime which is reached asymptotically, when the non linear saturation is observed, is characterized by the natural Strouhal number St n that can be estimated by a spectral analysis of the aerodynamic coefficients. On Fig.8, it can be observed that the drag force consists only of contributions from the even harmonics, and the lift force of contributions from the odd harmonics only. This well-known empirical fact was recently explained in Protas and Wesfreid 34 by considering the symmetry properties of the global modes known to exist in periodic wake flows. In table II, the natural Strouhal number and the time-averaged drag coefficient are compared to reference results available in the literature. The agreement with all the previous experimental and computational data is very good. Similarly (not shown in tableII), the time-averaged lift coefficient is seen to be in very good agreement with the results obtained previously. We conclude our presentation of the results obtained at Re = 200 by a discussion of the base flow (the unstable symmetric state). Recently, Protas and Wesfreid 25 argued that the mean drag C D consists of the contribution of two terms: the drag C base D of the basic flow (i.e., the unstable, steady, symmetric flow) which at a given Reynolds number remains fixed and the drag C 0 D of the mean flow correction which is due to the vortex shedding:

  0 and u m (x) matches all other non homogeneous boundary conditions, the POD basis functions φ i satisfy homogeneous boundary conditions on the whole domain. As a consequence, there is no contribution of the pressure term in the POD ROM. Due to the non reflecting boundary conditions (5) used in the outflow

Figure 14

 14 Figure 14 represents schematically the above algorithm.

  the spectrum of the excitation function γ e is represented. The frequencies vary continuously from St = 0.15 to St ≃ 0.65 and the spectrum presents a weak dominating mode for St ≃ 0.4. The Navier-Stokes equations are then solved with γ e for boundary conditions on the cylinder. During the course of the excitation 600 snapshots are taken uniformly over one period T e = 60 of excitation. These snapshots are used to form the temporal correlation matrix for the Fredholm equation (11).The POD eigenvalues for the uncontrolled flow (γ = 0) and controlled flow (γ = γ e ) are shown in Fig.17on a semi log-scale. For the uncontrolled flow, the set of eigenvalues fall-off rapidly, and hence a low number of POD modes is necessary to represent accurately any velocity field. Clearly for the controlled flow the spread of energy is much more uniform and many more degrees of freedom than for the uncontrolled flow are excited. As a consequence for a given number of modes kept in the POD expansion, the projection error of the snapshots on the POD basis functions is greater for the controlled flow than for the uncontrolled flow. timization horizon considered in Ref. 21 where an optimal control approach based on the Navier-Stokes equations is developed to control almost the same configuration of wake flow.

2 ;

 2 St f = 0.53), do not correspond to those of the mean drag minimization that can be obtained by an open loop control study (A min = 4.3; St f min = 0.74), see Fig. 31 in Appendix B. However the wide "valley" corresponding to the minimal mean drag is determined and 88% of the relative drag minimization that can be numerically found for this flow configuration under rotary control is nevertheless obtained. Now, in Fig. 27 we represent the comparison of the time evolution of the lift coefficient for the uncontrolled and optimally controlled flow obtained when ℓ 2 = 0. We observe that the action of control reduces considerably the amplitude of the lift coefficient (from 1.38 to 0.34). These behavior are synthesized in Fig. 28 where the polar curves (time evolution of the drag coefficient versus the lift coefficient) are represented for the uncontrolled and optimally controlled flow (ℓ 2 = 0). The limit cycles appearing in this figure are well defined because each aerodynamic coefficient oscillates with only one frequency. The power spectral density of the corresponding aerodynamic coefficients represented in Fig. 29 demonstrate that the controlled flow now oscillates at the frequency of the optimal control function (St f = 0.53). Finally, in Figs. 30(a), 30(b) and 30(c) we represent the vorticity fields of the uncontrolled flow, the optimally controlled flow (ℓ 2 = 0) and the base flow respectively. The significant vortex-shedding phenomenon observed in Fig. 30(a) has been substantially reduced when the control is applied and the flow has been quasi-symmetrized. The resulting flow approaches the symmetric state characteristic of the corresponding base flow as it can be awaited from the discussion in Sec. II B. Our results are qualitatively similar to the effects observed in Refs. 19,24 and confirm the arguments of Protas and Styczek 21 that the mean drag reduction is associated with control driving the mean flow toward the unstable state.

  Sec. VI B), the use of a POD ROM for state equation permits to consider an optimization horizon equal to six times the natural vortex shedding. The research of He et al. 19 shows a 30% drag reduction if one uses a sinusoidal rotating cylinder with the amplitude A = 3 and the forcing Strouhal number St f = 0.75. For the same flow configuration and only a slightly lower value of the Reynolds number (Re = 100), Homescu et al. 20 found that the optimal control parameters are A = 3.25 and St f = 1.13. These optimal amplitude and forcing Strouhal number differ considerably from the values found with our approach. As it was already suggested by Homescu et al., these difference appear to be due to the different formulations of the cost functionals used by the different authors. In Ref. 20, the main goal was the suppression of the Kármán vortex shedding and the cost functional was chosen to be of the flow tracking type (the "desired" flow correspond to a creeping flow at Re = 2), while in Ref. 19, He et al. aimed to reduce the drag and minimize a drag-related cost function (see their equation 9) and, while in Ref. 21, Protas and Styczek considered as cost functional directly the sum of the work done against the drag force and the work needed to control the flow.

FIG. 4 .

 4 FIG. 4. Configuration of controlled flow.

FIG. 10 .

 10 FIG. 10. Streamlines for the base flow at Re = 200.

  FIG. 1: Bergmann, Cordier and Brancher, Physics of Fluids

γFIG

  FIG.14:

(a)A = 3 ;FIG

 3 FIG. 32: Bergmann, Cordier and Brancher, Physics of Fluids

  ) are integrated in time, numerical instabilities arise after a few vortex shedding period (see Ref. 36 for example) and the model is no longer sufficiently accurate. This problem is similar to that of Large Eddy Simulation where we have to model the energy transfers between the Fourier modes lower than a given cutoff value that are simulated and those higher than this cutoff value that are not explicitly simulated. Recently, Karniadakis employed the same dissipative model called Spectral Vanishing Viscosity Model (SVVM) to formulate alternative LES approaches 37

)

  FollowingDeane et al. 36 and Noack et al. 14 , this term is found to be negligible for a wake flow configuration. Recently,Noack et al. 46 used a modal energy-flow analysis to elucidate the effect of the pressure term in a POD ROM of incompressible shear-flows. Essentially, they demonstrated that the effect of the pressure term is important for a mixing layer and small in a wake flow. In our approach, the pressure term (A3) is indeed omitted in the coefficients of the POD ROM (see coefficients A i in Sec. A 2). Here, this omission is numerically justified by the introduction of a time dependent eddy-viscosity for each POD mode (the method is described inRef. 32). With this dissipative model, an accurate low order modelling of a controlled wake flow is possible (see Figs.12 and 13) and neglecting the pressure term in the linear coefficients A i have no influence.

  Fahl presents in Ref. 45 an algorithm that implements the combination of POD based reduced order modelling and trust region methods, the TRPOD (Trust Region Proper Orthogonal De-

	composition).

TABLE II :

 II . 1. Discussion of parameter sampling in an optimization setting (figure taken from Gunzburger 51 ). ---path to optimizer using full system, initial values used by the optimizer and optimal values, • parameter values used for snapshot generation.

	List of figures
	Re	Authors	St n	C D
	200	Braza et al. 49	0.2000	1.4000
		Henderson 50	0.1971	1.3412
		He et al. 19	0.1978	1.3560
		Present study	0.1999	1.3900
	FIG. 2. Conceptual optimization algorithm using POD Reduced-Order Model.
	Bergmann, Cordier and Brancher, Physics of Fluids FIG. 3. Conceptual optimization algorithm using POD Reduced-Order Model (adaptive
	method).

FIG

Acknowledgments

The authors acknowledge M. Braza (Institut de Mécanique des Fluides de Toulouse) and D. Ruiz (ENSEEIHT) to kindly provide us with an original version of their Matlab Navier-Stokes solver. Laurent Cordier acknowledges the low-dimensional modelling and control team at the Technische Universität Berlin, in particular Bernd R. Noack and Ivanka Pelivan, for fruitful and stimulating discussions during the process of this work. Finally, we thank anonymous referees for their valuable suggestions.

Sec. IV. This quantity is represented in Fig. 18 for the uncontrolled (γ = 0) and controlled (γ = γ e ) flows. It is found that two modes are sufficient to represent 98% of the total kinetic energy in the uncontrolled case and that it is necessary to keep 40 modes (less than 7% of all the POD modes) to capture the same percentage of energy when the excitation γ e is applied. Therefore the state equation ( 22) used in the iterative optimization process of Sec. V correspond to N gal = 40.

Velocity contours for the first six POD eigenfunctions are plotted in Fig. 19. A striking feature is the pairing (mode 1 and 2, mode 4 and 5) of similar patterns, shifted spatially, a result of the convective nature of the flow. This behavior can be noticed in Fig. 17 where the corresponding eigenvalues occur in pairs of almost equal values. Clearly, mode 3 is localized around the cylinder and corresponds to the action of the control on the flow.

In first approximation, this boundary-layer mode is comparable to the control function u c determined previously (see Fig. 11 for comparison).

B. Results of the POD ROM based control

In this section and in the rest of this paper, we consider the results of the optimal control approach described in Sec. V for two couples of regularization parameters (see Eq. 24) i.e.

(ℓ 1 ; ℓ 2 ) = (1 ; 0) and (ℓ 1 ; ℓ 2 ) = (1 ; 1). When ℓ 2 = 0, it corresponds to the case where the cost of the control is neglected. Since, our main concern in the present study is not the energetic efficiency of the control procedure, only this case will be discussed in details.

Of course, as it can be expected and obtained (see Sec. VI D), when the control input is penalized, the energetic efficiency of the algorithm improves considerably and lower values of the amplitude A and the forcing Strouhal number St f are found (see hereafter). For higher values of the penalization parameter ℓ 2 , the optimization process does not converge any more. This divergence has certainly to be related to the choice of the parameters used to define the temporal excitation γ e . The robustness of the POD basis functions determined with this specific excitation is not sufficient to describe the controlled dynamics for rather small amplitudes of the control parameters.

Since the excitation function γ e is symmetric with respect to t = 30, the optimization horizon T o is restricted to T o = 30. This value corresponds to 5 times the maximum op-ported in Refs. 19,20,24,25,43. Therefore, if we want to compare our results with those of the literature, we need to determine the amplitude A and the forcing Strouhal number St f such that the optimized control function writes γ opt (t) = A sin(2πSt f t). A time average amplitude A ≃ 2.2 is easily determined. In Fig. 25 we show the spectrum of the control function γ opt obtained for ℓ 2 = 0: a main peak appears corresponding to St f ≃ 0.53 and two other peaks much less energetic are also clearly visible. According to He et al. 19 the contributions of these lower modes are negligible for the drag. Here, this point was checked by using this sinusoidal law as an open-loop control function. The relative reduction of the wake unsteadiness J varies from a value equal to 69% when no particular assumption is done on the variation of the control function to a value equal to 52% when the control function is considered sinusoidal. The relative reduction is still consequent. Therefore the two other peaks are neglected in the following and the optimal control γ opt is supposed to be sinusoidal with an amplitude A = 2.2 and a Strouhal number St f = 0.53 when ℓ 2 = 0.

When ℓ 2 = 1, exactly the same numerical procedure was adopted. As it can be intuitively expected, we obtain lower values for the amplitudes of the control parameters, respectively A = 0.75 and St f = 0.33.

C. Drag reduction for the Navier-Stokes model

By definition of the optimization problem (25), the control function γ opt is optimal for the POD ROM and not necessary for the Navier-Stokes model. More exactly, with this particular approach, there is no mathematical assurance that the iterates produced by the optimization algorithm will converge to a local optimizer for the high-fidelity original problem (certainly the trust region idea from nonlinear programming 44 would be fruitful to circumvent this difficulty). Of course, it was found in Sec. VI B that the optimized control function γ opt reduced considerably the wake unsteadiness, but the initial objective of this paper is the drag minimization. Therefore it is necessary to solve the Navier-Stokes equations with a rotary control defined by γ opt (t) to determine the effect of this control function on the drag coefficient. When the POD ROM ( 22) is derived from the velocity expansion (18) with a standard Galerkin projection on the Navier-Stokes system (1), it appears a contribution from the pressure term:

where

Since the POD basis functions φ i can be represented as linear combinations of instantaneous velocity fields (see Eq. 14), they inherit all the properties of the snapshots that can be written as linear and homogeneous equations. So, if the velocity fields included in the input data are solenoidal, then divergence-free POD basis functions are obtained (∇ • φ i = 0) and (A1) becomes:

When the control function method is employed, POD basis functions with homogeneous boundary conditions (φ i = 0) are determined on all the boundaries where Dirichlet conditions were considered in the numerical simulation (see Sec. IV A). Therefore, there is no pressure contribution from the inflow and cylinder boundaries. Moreover, due to the boundary conditions prescribed on the top and bottom walls (see Eq. 4), the POD functions verify the condition φ i • n = 0 and the contributions from the top and bottom boundaries are then A series of simulations with different amplitude A varying from 0 to 6.5 by step of 0.5 and different forcing Strouhal number St f varying from 0 to 1. by step of 0.1 was performed.

natural shedding frequency St n . For every forcing frequency our simulations are performed for a sufficient long time (T S = 130) to assure that the saturated state has been reached.

All simulations have been done with the same time step, here equal to 1.5 10 -2 . In Fig. 31, we visualize the contours of the mean temporal drag estimated over the last 30 units of time as a function of A and St f . Numerically, the mean drag is minimized for an optimal pair (A min , St f min ) = (4.3, 0.74). The corresponding minimum value is 0.99. Finally, to illustrate the different flow patterns that can be obtained for the forced flow, we represent in Fig. 32 the vorticity fields corresponding to a fixed value equal to 3 of the amplitude A and different forcing Strouhal numbers. to control the laminar wake flow with an optimal control approach. The present study is included for comparison. 'Unknown' means that the value was not found in the article, an estimation is then given when it is possible.

List of tables