Model-based Stabilization of Vortex Shedding with CFD Veri cation

Milan Milovanovic, Marek Gayer, Joris Michielsen and Ole Morten Aamo

Abstract— We establish the link between backstepping con-

trollers previously designed for the Ginzburg-Landau model of

vortex shedding and the actual ow dynamics, represented by

a Navier-Stokes solver, and demonstrate, for the rst time in

CFD simulations, the effectiveness of this class of controllers in

attenuating vortex shedding.

I. I

NTRODUCTION

In ows past submerged obstacles, the phenomenon of

vortex shedding occurs provided the Reynolds number is

suf ciently large. A popular model ow for studying vortex

suppression by means of open-loop or feedback control, is

the ow past a two-dimensional circular cylinder, as sketched

in Fig. 1 where the the so-called von Kármán vortex street

is visualized using passive tracer particles. For Reynolds

numbers slightly larger than the critical value for onset

of vortex shedding (which is approximately Re

c

= 47),

several authors have successfully suppressed vortex shedding

in numerical simulations using various simple feedback

control con gurations. In [1], a pair of suction/blowing slots

positioned on the cylinder wall were used for actuation, and

shedding was suppressed for Re = 60, using proportional

feedback from a single velocity measurement taken some

distance downstream of the cylinder. For Re = 80, vortex

shedding was reduced, but not completely suppressed. In [2],

the same actuation con guration was tried using feedback

from a pair of pressure sensors located on the cylinder wall

for Re = 60. This attempt was unsuccessful, but by adding

a third actuation slot shedding was reduced considerably,

even at Re = 80. Although some success in controlling

vortex shedding has been achieved in numerical simulations,

rigorous control designs are scarce due to the complexity of

designing controllers based on the Navier-Stokes equation.

A much simpler model, the Ginzburg-Landau (GL) equation

with appropriate coef cients, has been found to model well

the dynamics of vortex shedding near the critical value of

the Reynolds number [3]. In [4], it was shown numerically

that the GL model for Reynolds numbers close to Re

c

can be stabilized using proportional feedback from a single

measurement downstream of the cylinder, to local forcing

at the location of the cylinder. Systematic control designs

for nite dimensional approximations of the model from [4]

were developed in [5], [6] for the linearized model, and in [7]

for the nonlinear model. The designs were in principle valid

for any Reynolds number. In [8] and [9], boundary control

This work was supported by the BVV program at NTNU (Program on

Computational Science and Visualization). The second author acknowledges

support from the ERCIM exchange program.

M. Milovanovic, M. Gayer, and O.M. Aamo are with the Department of

Engineering Cybernetics, NTNU. J. Michielsen is with TU Eindhoven.

Fig. 1.

Vortex shedding from a cylinder visualized by passive tracer

particles

laws based on the non-discretized, linearized GL equation

were rigorously derived for the state feedback and boundary

output feedback cases, respectively. They were shown to

effectively attenuate vortex shedding for the nominal case

with GL as the system plant as well as the design model.

The link between the GL-based controller and the physical

system, represented by a Navier-Stokes (NS) solver, has not

yet been established. Closing this gap is the main topic

of this paper, which includes tailoring a CFD code for

ow control simulations, tting the GL model to CFD data,

and designing the input-output links between the GL-based

controller and the NS-based simulator. Putting all pieces

together, we demonstrate in CFD simulations the successful

stabilization of vortex shedding by means of the GL-based

controller.

II. T

HE

G

INZBURG

-L

ANDAU

M

ODEL

The dynamics of the cylinder wake is governed by the

Navier-Stokes equation. In [4], however, a simpli ed model

was suggested in terms of the Ginzburg-Landau equation

@A

@t

=

a

1

@

2

A

@x

2

+ a

2

(x)

@A

@x

+ a

3

(x) A

+a

4

jAj

2

A + (x

1) u;

(1)

where A is a complex-valued function of one spatial variable,

x 2 R, and time, t 2 R

+

. The boundary conditions are

A ( 1;t) = 0. The control input, denoted u, is in the form

of point actuation at the location of the cylinder, and the

coef cients a

i

; i = 1; ::; 4; were tted to data from laboratory

experiments in [4].

denotes the Dirac distribution. A (x; t)

may represent any physical variable (velocities (u; v) or

pressure p), or derivations thereof, along the centerline of

the 2D cylinder ow. The choice will have an impact on the

performance of the Ginzburg-Landau model, and associating

A

with the transverse

uctuating velocity v (x; y = 0; t)

seems to be a particularly good choice [10]. As pointed out in

[6], the model is derived for Reynolds numbers close to the

critical Reynolds number for onset of vortex shedding, but

has been shown to remain accurate far outside this vicinity

for a wide variety of ows.

In [8], [9] a simpli cation of (1) was considered, ob-

tained by linearizing around the zero solution, discarding

the upstream subsystem by replacing the local forcing at

x = 0 with boundary input at this location, and truncating

the downstream subsystem at some x

d

2 (0; 1). The

truncation of the system was justi ed by noting that the

upstream subsystem is approximately uniform ow, whereas

the downstream subsystem can be approximated to any

desired level of accuracy by selecting x

d

suf ciently far from

the cylinder.

1

The resulting system is given by

@A

@t

= a

1

@

2

A

@x

2

+ a

2

(x)

@A

@x

+ a

3

(x) A

(2)

for x 2 (0;x

d

); with boundary conditions

A (0; t) = u(t) and A (x

d

; t) = 0;

(3)

where A : [0; x

d

]

R

+

! C, a

2

2 C

2

([0; x

d

]; C) ;

a

3

2 C

1

([0; x

d

]; C), a

1

2 C, and u : R

+

! C is the

control input. a

1

is assumed to have strictly positive real

part. In [8], stabilizing state feedback boundary control laws

for system (2)–(3) were derived based on the backstepping

methodology [11], [12], [13], [14]. The control laws made

use of distributed measurements in a nite region down-

stream of the cylinder. In practice, the physical actuation

could be for instance micro/synthetic jet actuators, distributed

on the cylinder surface, and a possible choice for the physical

sensing could be pressure sensors distributed on the cylinder

surface.

III. S

TATE

F

EEDBACK

C

ONTROL

D

ESIGN

The state feedback controller designed in [8] takes the

form

u(t)

=

Z

x

d

0

1

x

d

[k

1

(

x

d

x

x

d

)

ik

c;1

(

x

d

x

x

d

)]A(x; t)

exp

1

2a

1

Z

x

0

a

2

( )d

!dx

(4)

where i denotes the imaginary unit, and k

1

and k

c;1

are the controller kernels which can be computed as the

solution of a certain hyperbolic partial differential equa-

tion. We repeat some details from [8] for clarity and

completeness. It is convenient to rewrite the GL equa-

tion to obtain two coupled partial differential equations

in real variables and coef cients by de ning

(x; t) =

<(B(x;t)) = B (x;t) + B (x;t) =2, and (x;t) =

=(B(x;t)) = B (x;t) B (x;t) =(2i), where x =

(x

d

x) =x

d

; B (x; t) = A (x; t) exp

1

2a

1

R

x

0

a

2

( ) d

;

and

denotes complex conjugation. Equation (2) becomes

t

=

a

R xx

+ b

R

(x)

a

I xx

b

I

(x) ;

t

=

a

I xx

+ b

I

(x) + a

R xx

+ b

R

(x) ;

(5)

1

This claim is postulated from the observation that the local damping

effect in (1) increases with increasing distance from the cylinder, which

follows from the coef cients reported in [4].

for x 2 (0;1), with boundary conditions

(0; t)

=

0;

(0; t) = 0; and

(1; t)

=

u

R

(t) ;

(1; t) = u

I

(t) ;

(6)

where a

R

, <(a

1

)=x

2

d

; a

I

, =(a

1

)=x

2

d

;

and

b

R

(x)

, < a

3

(x)

1

2

a

0

2

(x)

1

4a

1

a

2

2

(x) ;

b

I

(x)

, = a

3

(x)

1

2

a

0

2

(x)

1

4a

1

a

2

2

(x) :

(7)

The state feedback stabilization problem is solved by search-

ing for a coordinate transformation in the form

(x; t) =

(x; t)

Z

x

0

[k (x; y) (y; t) + k

c

(x; y) (y; t)] dy;

(8)

(x; t) = (x; t)

Z

x

0

[ k

c

(x; y) (y; t) + k (x; y) (y; t)] dy;

(9)

transforming system (5)–(6) into

t

=

a

R xx

+ f

R

(x)

a

I xx

f

I

(x) ;

t

=

a

I xx

+ f

I

(x) + a

R xx

+ f

R

(x) ;

(10)

for x 2 (0;1), with boundary conditions

(0; t) = (0; t) =

(1; t) = (1; t) = 0:

(11)

By the choice of f

R

and f

I

, system (10)–(11) can be

given any desired level of stability. The corresponding stable

behaviour for the original system is ensured by the control

input

u

R

(t) =

Z

1

0

[k

1

(y) (y; t) + k

c;1

(y) (y; t)] dy;

(12)

u

I

(t) =

Z

1

0

[ k

c;1

(y) (y; t) + k

1

(y) (y; t)] dy;

(13)

where k

1

(y) = k (1; y) ; k

c;1

(y) = k

c

(1; y). The skew-

symmetric form of (12)–(13) is postulated from the skew-

symmetric form of (5). The following result was proven in

[8]

2

.

Theorem 1:

i. The pair of kernels, k (x; y) and k

c

(x; y), satisfy the

partial differential equation

k

xx

=

k

yy

+ (x; y)k +

c

(x; y)k

c

;

k

c;xx

=

k

c;yy

c

(x; y)k + (x; y)k

c

;

(14)

for (x; y) 2 T = fx;y : 0 < y < x < 1g, with boundary

conditions

k(x; x)

=

1

2

R

x

0

( ; )d ; k (x; 0) = 0;

k

c

(x; x)

=

1

2

R

x

0

c

( ; )d ; k

c

(x; 0) = 0;

(15)

where

(x; y) =

a

R

(b

R

(y)

f

R

(x)) + a

I

(b

I

(y)

f

I

(x))

a

2

R

+ a

2

I

;

(16)

2

In the statement of Theorem 1 and in the text in Section 8 of [8],

L

1

(0; 1)

should be replaced by H

3

(0; 1)

. Furthermore, the initial condi-

tions must be compatible with the boundary conditions.

c

(x; y) =

a

R

(b

I

(y)

f

I

(x))

a

I

(b

R

(y)

f

R

(x))

a

2

R

+ a

2

I

:

(17)

The equation (14) with boundary conditions (15) has a

unique C

2

(T ) solution, given by

k (x; y)

=

P

1

n=0

G

n

(x + y; x

y) ;

k

c

(x; y)

=

P

1

n=0

G

c;n

(x + y; x

y) ;

(18)

where

G

0

( ; )

=

1

4

Z b( ;0)d ; (19)

G

c;0

( ; )

=

1

4

Z b

c

( ; 0) d ;

(20)

G

n+1

( ; ) =

1

4

Z Z

0

b ( ; s) G

n

( ; s) dsd

+

1

4

Z Z

0

b

c

( ; s) G

c;n

( ; s) dsd ;

(21)

G

c;n+1

( ; ) =

1

4

Z Z

0

b

c

( ; s) G

n

( ; s) dsd

+

1

4

Z Z

0

b ( ; s) G

c;n

( ; s) dsd ;

(22)

and

b ( ; ) =

+

2

;

2

;

(23)

b

c

( ; ) =

c

+

2

;

2

:

(24)

ii. The inverse of (8)–(9) exists and is in the form

(x; t) =

(x; t)

Z

x

0

[l (x; y) (y; t) + l

c

(x; y) (y; t)] dy;

(25)

(x; t) = (x; t)

Z

x

0

[ l

c

(x; y) (y; t) + l (x; y) (y; t)] dy;

(26)

where l and l

c

are C

2

(T ) functions. l and l

c

can be expressed

similarly to k and k

c

in (18)–(22), but we omit their explicit

de nition due to page limitations.

iii. Suppose c > 0, and select f

R

and f

I

such that

sup

x

2[0;1]

f

R

(x) +

1

2

jf

0

I

(x)j

c:

(27)

Then for any initial data (

0

;

0

) 2 H

3

(0; 1), compatible

with the boundary conditions, the system (5)–(6) in closed

loop with the control law (12)–(13) has a unique classical

solution ( ; ) 2 C

2;1

((0; 1)

(0; 1)) and is exponentially

stable at the origin in the L

2

(0; 1) and H

1

(0; 1) norms.

In [8], it was shown that a particular choice of f

R

and

f

I

, that depend on x

d

in a speci c way, results in state

feedback kernel functions that are invariant of x

d

and vanish

Fig. 2. Simulation setup.

in [x

s

; x

d

], where x

s

is a constant that can be deduced from

the coef cients of (2). This implies that if the domain is

truncated at some x

s

x

d

for the purpose of computing

the feedback kernel functions, the resulting state feedback

will stabilize the plant evolving on the semi-in nite do-

main (0; 1). This is achieved by avoiding the complete

cancellation of the terms involving b

R

and b

I

in (5) by

using a target system (10) that contains the natural damping

that exists in the plant downstream of x

s

. It ensures that

only cancellation/domination of the source of instability is

performed in the design, and is similar to common practice in

design of nite dimensional backstepping controllers, where

one seeks to leave unaltered terms that add to the stability

while cancelling terms that don't. The result is less complex-

ity, and better robustness properties. An interesting property

of our design is that it requires the solution of a linear

hyperbolic PDE, which is an advantage when compared to

other methods, such as LQG requiring the solution of a

Riccati equation, which is quadratic. In fact, for a plant

much simpler than the linearized Ginzburg-Landau model,

solving the hyperbolic PDE is reported in [11] to take an

order of magnitude less computational time than solving the

Riccati equation. Although the series (18) happen to converge

rapidly, they can be computed by a more ef cient numerical

scheme that was developed in [11]. Further details can be

found in [8].

IV. S

IMULATION

S

ETUP

The simulation setup is shown schematically in Fig. 2.

It contains of ine computations performed in MATLAB that

produce GL parameters for a given Reynolds number and the

corresponding controller kernels. The online computations,

or the closed-loop simulations, are performed with the CFD

code VISTA. We next present details of the simulation setup.

-0.5

0

0.5

1

1.5

2

-1

-0.5

0

0.5

1

Fig. 3. Curves along which transverse velocity is measured.

A. GL/CFD Connection

While A in principle may represent any physical variable

in the ow, the transverse uctuation velocity along the ow

centerline is a particularly good choice. For the interior ow,

we de ne two curves originating at the actuation slots and

converging at the ow centerline as they stretch into the 2D

ow domain (see Fig. 3). Transverse velocity is measured

along these curves, denoted v

l

(x; t) and v

u

(x; t) in Fig. 3,

and related to A as

Re(A) := v(x; t) =

1

2

(v

l

(x; t) + v

u

x; t)) :

(28)

The imaginary part of A is reconstructed based on modes of

the GL equation having the form

A(x; t) =

(x)exp(i x

i!t):

(29)

We approximate Im(A) by shifting Re(A) in space corre-

sponding roughly to

x = T U=4, where T is the shedding

period and U is the streamwise velocity. Having obtained A

as described, the control input is computed using Equation

(4). With the real part of A representing transverse velocity,

the control actuation which is a Dirichlet boundary condition

at the actuation end in the GL model, can be realized in the

CFD code by jets on the cylinder surface. We do this by

employing two slots con gured as shown in Fig. 4, close to

the separation layers for effectiveness [15]. When one slot

applies suction the other applies an equal amount of blowing,

thereby maintaining the mass balance in the system. The

actuation is thus suction and blowing whose size is given by

the real part of u computed from Equation (4).

B. Obtaining GL parameters

To estimate the parameters for the Ginzburg-Landau equa-

tion, an iterative search is carried out where the difference

between the model solution and a data set is minimized.

The data set is extracted from CFD simulations, where

the transition between laminar

ow and fully developed

vortex shedding is made by changing the Reynolds number

from subcritical (R

e

= 40) to supercritical (R

e

= 60).

The transverse velocity component, v(x; t), is measured as

described in the previous section along the curves shown in

*x*

*y*

Fig. 4. Suction/blowing actuation slots.

Fig. 5. CFD data set.

Fig. 3. The resulting data set is shown in Fig. 5, where the

transition into vortex shedding is clearly visible.

By studying CFD data and simulations of the GL equation

(3), some apriori assumptions were made for the parameters

a

1

, a

2

(x) and a

3

(x) before posing the optimization problem.

The diffusion contribution is assumed to be constant, a

1

=

1=R

e

, because it corresponds to the dynamic viscosity of

the uid. a

2

(x) is assumed to be real-valued, while a

3

(x) =

a

3R

(x)+ia

3I

(x) is complex-valued with components a

3R

(x)

and a

3I

(x) real-valued. The functions a

2

(x), a

3R

(x) and

a

3I

(x) are approximated as polynomials in the optimization

problem.

To obtain an accurate estimate for the parameters, a

Fig. 6. Estiamated parameters a

2;

a

3R;

a

3I

ordered from top to bottom

Fig. 7. The residue between real CFD data and GL simulations.

0

0.5

1

1.5

2

-50

-40

-30

-20

-10

0

10

20

30

x

Fig. 8. Controller kernels, k

1

(solid line) and k

c;1

(dashed line).

nite difference scheme for the GL equation is derived

based on the knowledge that the system is dominated by

advection. To avoid numerical diffusion, which leads to

severe damping in the advection term, the Lax-Wendroff

scheme is used. This scheme is second order accurate in

space and time. To keep this accuracy a central difference

scheme for the spatial second order derivative in the diffusion

term is combined with the Crank-Nicholson scheme as time

integration method for both the diffusion and advection

terms. Denoting the solution of the difference scheme on the

spatial grid fx

1

; x

2

; :::; x

m

g as A

R

= [A

R

(x

1

) ::: A

R

(x

m

)]

T

and A

I

= [A

I

(x

1

) ::: A

I

(x

m

)]

T

(corresponding to the real

and imaginary parts), the optimization method is based on

the least squares formulation

c = arg min

c

f (c) = arg min

c

1

N

N

X

n=1

1

2

r

n

(c)

T

r

n

(c)

where r

n

(c) = A

n

R

(c)

v

n

and n speci es the discrete

time step, c is a vector with unknown parameters to be

determined, and v

n

is the CFD data at time step n. This

unconstrained nonlinear optimization problem is solved with

the Levenberg-Marquardt method.

For the data set shown in Fig. 5, the initial parameter guess

and the converged parameters are presented in Fig. 6. The

error when comparing the data set with simulations with the

GL equation using these paramters is shown in Fig. 7. The

parameters are effetively estimated.

With parameters of the Ginzburg-Landau equation in hand,

the controller kernels, k

1

(x) and k

c;1

(x), are computed from

(18)–(24). They are shown in Fig. 8.

Fig. 9. Computational grid with zoom around cylinder.

Fig. 10. Vorticity map at t = 0: Fully developed vortex shedding.

C. CFD system

As a basis for our CFD simulations we use VISTA,

which is a general purpose CFD code developed by SINTEF

Applied Mathematics. At the core of VISTA is a Navier-

Stokes solver written in C++ using the Diffpack library for

nite element modelling. We have tailored VISTA by adding

C++ objects that realize actuation and sensing devices, as

well as controller computations. This module takes as input

the controller kernels that are computed of ine in MATLAB.

General purpose grid generators can be used to produce

computational grids for VISTA. Fig. 9 shows the grid used in

our simulations, which was veri ed in [16] to be suf ciently

ne for the Reynolds numbers that we use in this paper.

V. S

IMULATION

R

ESULTS

The intial condition for our simulation is fully developed

vortex shedding at R

e

= 60. Fig. 10 shows a vorticity

map of the ow at t = 0. The von Kármán vortex street

is clearly visible. Vortex shedding subjects the bluff body,

in this case the cylinder, to periodic forcing, and would in

any physical system induce potentially destructive vibrations.

The transverse force is shown in Fig. 11, for our controlled

simulation. The controller effectively attenuates the forcing.

Fig. 12 shows the input to our controller, Re(A), as a

function of time, while Fig. 13 and Fig. 14 show vorticity

maps at t = 20s and t = 100s. Clearly, vortex shedding is

completely removed in Fig. 14. The control actuation that

achieves this is shown in Fig. 15.

VI. C

ONCLUSION

In this paper, we have established the link between

controllers previously designed for the Ginzburg-Landau

0

20

40

60

80

100

-0.2

-0.1

0

0.1

0.2

Fig. 11. Lift coef cient.

Fig. 12. Input to the controller, A(x; t):

model of vortex shedding and the actual ow dynamics. The

latter has been represented by a CFD package that solves

the Navier-Stokes system for the 2D cylinder ow. Vortex

shedding is successfully removed by the control at Reynolds

number R

e

= 60. The control design is systematic and

can potentially be repeated for higher Reynolds numbers,

although it is likely to be limited by the validity of the

Ginzburg-Landau model at high Reynolds number. Ongoing

work is focused at exploring this issue, while future work will

involve CFD veri cation of the output feedback problem,

which was solved for the GL equation in [9].

R

EFERENCES

[1] D. Park, D. Ladd, and E. Hendricks, “Feedback control of von kármán

vortex shedding behind a circular cylinder at low reynoylds numbers,”

Physics of uids, vol. 7, no. 7, pp. 2390–2405, 1994.

[2] M. Gunzburger and H. Lee, “Feedback control of karman vortex

shedding,” Transactions of the ASME, vol. 63, pp. 828–835, 1996.

Fig. 13. Vorticity map at t = 20s.

Fig. 14. Vorticity map at t = 100: Completely stabilized ow.

Fig. 15. The control input.

[3] P. Huerre and P. Monkewitz, “Local and global instabilities in spatially

developing ows,” Annu. Rev. Fluid Mech., vol. 22, pp. 473–537, 1990.

[4] K. Roussopoulos and P. Monkewitz, “Nonlinear modelling of vortex

shedding control in cylinder wakes,” Physica D, vol. 97, pp. 264–273,

1996.

[5] E. Lauga and T. Bewley, “H in nity control of linear global insta-

bility in models of non-parallel wakes,” Proceedings of the Second

International Symposium on Turbulence and Shear Flow Phenomena,

vol. Stockholm, Sweden, 2001.

[6] E. Lauga and T. Bewley, “Performance of a linear robust control

strategy on a nonlinear model of spatially-developing ows,” J. Fluid

Mech, vol. 512, pp. 343–374, 2004.

[7] O. M. Aamo and M. Krstic, “Global stabilization of a nonlinear

ginzburg-landau model of vortex shedding,” Eur. J. Control, vol. 10,

no. 2, 2004.

[8] O. M. Aamo, A. Smyshlyaev, and M. Krstic, “Boundary control of

the linearized ginzburg-landau model of vortex shedding,” SIAM J.

Control Optim, vol. 43, no. 6, pp. 1953–1971, 2005.

[9] O. M. Aamo, A. Smyshlyaev, M. Krstic, and B. Foss, “Output

feedback boundary control of a ginzburg-landau model of vortex

shedding,” IEEE Transactions on Automatic Control, vol. 52, no. 4,

pp. 742–748, 2007.

[10] P. Monkewitz, private communication.

[11] A. Smyshlyaev and M. Krstic, “Closed form boundary state feedbacks

for a class of 1d partial integro-differential equations,” IEEE Transac-

tions on Automatic Control, vol. 49, pp. 2185–2202, 2004.

[12] W. J. Liu, “Boundary feedback stabilizaton of an unstable heat

equation,” SIAM J. Control Optim, vol. 42, pp. 1033–1043, 2003.

[13] M. Krstic and A. Smyshlyaev, Boundary Control of PDEs: A Course

on Backstepping Designs. SIAM, 2008.

[14] M. Krstic, I. Kanellakopoulos, and P. Kokotovic, Nonlinear and

Adaptive Control Design. Wiley, 1995.

[15] P. Catalano, M. Wang, G. Iaccarino, I. F. Sbalzarini, and P. Koumout-

sakos, “Optimization of cylinder ow control via zero net mass ux

actuators,” in Proceedings of the CTR summer program, Center for

Turbulence Research, Stanford University, 2002.

[16] P. Hoeijmakers, “Implementation of an extremum seeking controller

for vortex shedding attenuation in a 2d cfd code.” TU Eindhoven,

2008.