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.
