Optimal Dirichlet boundary control
problems of high-lift configurations with
control and integral state constraints
vorgelegt von
Diplom-Wirtschaftsmathematiker
Dipl.-Math.oec. Christian John
aus Berlin
Von der Fakultät II - Mathematik und Naturwissenschaften
der Technischen Universität Berlin
zur Erlangung des akademischen Grades
Doktor der Naturwissenschaften
Dr. rer. nat.
genehmigte Dissertation
Promotionsausschuss:
Vorsitzender: Prof. Dr. John M. Sullivan
Gutachter: Prof. Dr. Fredi Tröltzsch
Gutachter: Prof. Dr. Arnd Rösch
Tag der wissenschaftlichen Aussprache: 13.07.2011
Berlin 2011
D 83
ii
Abstract
This thesis investigates optimal control problems related to Navier-Stokes
equations. We investigate two control problems related to the aerodynamic
optimization of flows around airfoils in high-lift configurations.
The first issue is the steady state maximization of lift subject to restric-
tions on the drag. This leads to a Dirichlet boundary control problem for
the stationary Navier-Stokes equations with constrained control functions be-
longing to L2under an integral state constraint. The control space L2makes
it necessary to deal with very weak solutions of the Navier-Stokes equations
and because of the low regularity of control and state, we reformulate the cost
functional and the integral state constraint. We derive first-order necessary
and second-order sufficient optimality conditions and treat the problem nu-
merically by direct solution of the associated nonsmooth optimality system
and additionally by an SQP-method, which convergence we proved.
The second part is based on a k-ω-Wilcox98 turbulence model, describ-
ing the nonstationary behavior of the fluid closer to the reality. To deal with
the curse of dimension, we discuss a reduced-order model (ROM) by adapting
a small system of ODEs to solutions computed with the full model. Based
on this ROM, we investigate an optimal control problem theoretically and
numerically.
Acknowledgment
I think it is impossible to write a doctoral dissertation without some kind of
support. I am grateful to the DFG (SFB 557 ’Control of complex turbulent
shear flows’) supporting me in the first three years financially and to the
FAZIT-Stiftung for their financial support in the last 9 months.
My first and biggest gratitude goes to my supervisor Prof. Dr. Fredi
Tröltzsch for the interesting topic, his helpful advice and his personality in
general. I thank Dr. Daniel Wachsmuth, who worked with me together for
about 2 years when i started my project at the TU Berlin, for introducing me
into the topic of optimal control of Navier-Stokes equations, his interest in my
work and many helpful inspirations. Many thanks go to my colleagues in the
research group ’Optimization on PDEs’ at the TU Berlin, especially Kristof
Altmann for introducing me into COMSOL Multiphysics, and to Prof. Dr.
Arnd Rösch for his willingness to review this thesis.
iii
Furthermore, my gratitude to B.R. Noack and M. Schlegel for the good
and successful cooperation. They introduced me into the topic of proper or-
thogonal decomposition and reduced-order modeling and supported me with
both words and deeds. I also want to thank M. Luchtenburg for explaining
me his reduced-order model, M. Nestler for his helpful assistance and the
group of Prof. Thiele, especially B. Günther and A. Carnarius, for providing
me with simulations of the URANS system.
Finally, i thank my friends and my family.
iv
Contents
1 Introduction 1
2 The steady-state Navier-Stokes equation 7
2.1 Functionspaces.......................... 7
2.2 The Stokes equations . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 Weakformulation......................... 13
2.4 Very weak formulation . . . . . . . . . . . . . . . . . . . . . . 15
2.4.1 More regular solutions . . . . . . . . . . . . . . . . . . 17
2.4.2 Regularity assumption . . . . . . . . . . . . . . . . . . 19
3 The optimal control problem 21
3.1 Reformulation of the boundary integrals . . . . . . . . . . . . 21
3.2 Further reformulation of the boundary integrals . . . . . . . . 22
3.3 The optimal control problem . . . . . . . . . . . . . . . . . . . 24
3.3.1 Existence of solutions . . . . . . . . . . . . . . . . . . . 24
4 Optimality conditions 29
4.1 First order necessary optimality conditions . . . . . . . . . . . 29
4.2 Second-order sufficient optimality condition . . . . . . . . . . 33
4.3 Finite-dimensional control set . . . . . . . . . . . . . . . . . . 39
4.4 Second-order sufficient optimality conditions for the finite-
dimensionalcase ......................... 44
5 Numerical investigations 47
5.1 One-shotapproach ........................ 47
5.1.1 Numerical results . . . . . . . . . . . . . . . . . . . . . 49
5.2 SQP-method............................ 57
5.2.1 The SQP-method for our problem . . . . . . . . . . . . 58
5.2.2 Gradient-projection method . . . . . . . . . . . . . . . 65
5.2.3 Example.......................... 66
v
CONTENTS
6 Convergence of the SQP-method 69
6.1 Generalized equations . . . . . . . . . . . . . . . . . . . . . . . 71
6.2 Perturbed optimization problem . . . . . . . . . . . . . . . . . 76
6.3 Amodifiedproblem........................ 79
6.3.1 Existence of a solution . . . . . . . . . . . . . . . . . . 79
6.3.2 Lipschitz stability . . . . . . . . . . . . . . . . . . . . . 80
6.4 Strong regularity of the original perturbed problem . . . . . . 85
7 The nonstationary case 89
7.1 Modelreduction.......................... 91
7.2 Proper orthogonal decomposition POD . . . . . . . . . . . . . 92
8 Reduced-order model (ROM) 103
8.1 A generalized model . . . . . . . . . . . . . . . . . . . . . . . 106
8.1.1 Mean-field theory . . . . . . . . . . . . . . . . . . . . . 106
8.1.2 Mean-field Galerkin model . . . . . . . . . . . . . . . . 110
8.2 Modifications on the reduced-order model . . . . . . . . . . . 117
8.2.1 Filtering of the POD modes and coefficients . . . . . . 117
8.2.2 Parameter calibration . . . . . . . . . . . . . . . . . . . 124
8.3 Computation of lift . . . . . . . . . . . . . . . . . . . . . . . . 129
8.4 Numerical investigation . . . . . . . . . . . . . . . . . . . . . . 130
9 The optimal control problem 135
9.1 First-order necessary optimality conditions . . . . . . . . . . . 136
9.2 Numerical investigation . . . . . . . . . . . . . . . . . . . . . . 139
10 Conclusion 145
11 Zusammenfassung 147
vi
Chapter 1
Introduction
In this thesis, we study optimal control problems related to Navier-Stokes
equations, describing the motion of fluid. We investigate minimizations of
functionals subject to state equations. The objective functionals depend on
the velocity field u, the pressure pand the control function g.
Our main concern is maximization the lift of an airplane, while drag
remains beyond a given threshold. Therefore, we consider an objective func-
tional J(u, p, g)characterizing the lift, the Navier-Stokes equations as state
equations and a constraint on the drag.
In given literature, there are two different approaches to get influence on
the flow around a body. The first one is the possibility of passive control.
There are several possibilities of passive control, e.g. passive blowing,
roughness and shaping. Passive noise control devices include shields of rigid
and compliant walls, mufflers, silencers, resonators and absorbent materials,
see [43] for more details. The idea behind most of them is to reduce vortices
and make the airstream around the wing smoother.
The second ansatz is active flow control, which was investigated in partic-
ular by the SFB 557 ’Control of complex turbulent shear flows’. Here, little
slits are installed on a part of the wing, where suction and blowing of air is
possible to reduce vortices.
Generally, flow control is a research field gaining a lot of interest in both
academic research and industry. It is researched by engineers (experimen-
tal and computation fluid dynamics), mathematicians (control theory and
optimization) and physicists.
In this work, the following optimal flow control problem is considered: ac-
tive control of the flow of a fluid around an aircraft by means of suction and
blowing on the wing to influence the resulting lift and drag. The associated
background of applications in fluid mechanics, active separation control, was
the subject of various papers written from an engineering point of view and
1
CHAPTER 1. INTRODUCTION
has been proven to be effective in experiments as well as simulations. We
only mention [17, 19, 87, 88, 89, 112], whose considerations are close to our
setting, see Chapter 7 to 9.
The first part of this thesis deals with the steady-state problem. Here, we
assume a low Reynolds number so that we avoid the discussion of turbulence.
Furthermore, we consider a simplified control model, which is composed of
the cost functional, the steady-state Navier-Stokes equations, and constraints
on the control function as well as the state, for a mathematical investigation.
First, the steady-state Navier-Stokes equations, describing the motion of the
fluid around the wing, are investigated and we clarify the following questions.
1. What is the best (suitable) definition of a solution of the state equation
for the formulated problem?
2. What are the requirements such that a (unique) solution for the state
equations exists?
3. What preliminary results can be found in the existing literature?
4. What regularity assumptions are needed?
5. What are the requirements such that a solution for the stated opti-
mization problem exists?
6. How are the necessary and sufficient optimality conditions formulated?
7. What is an appropriate numerical optimization method and does it
converge?
We will characterize optimality of control strategies for our setting by
necessary and sufficient optimality conditions.
Let us now describe the setting of our optimization problem in detail.
Here, Ωis an open bounded domain of Rn,n= 2,3with boundary Γ, which
is assumed to be sufficiently smooth, more details later. The velocity field of
the fluid is denoted by uand the pressure by p. The control is a boundary
velocity field denoted by gand the viscosity parameter ν= 1/Re is a positive
number. Let us denote the surface measure by ds(x)or short ds. The
term ∇denotes the gradient and ∆the Laplace operator, which is applied
componentwise. The resulting force of the fluid on the wing embedded in the
fluid in direction ~e is given as the boundary integral
F~e =ZΓwν∂y
∂nw−pnw·~e ds, (1.0.1)
2
where Γwis the boundary of the wing with its outer normal nw. Since nw
points into the fluid, the normal nwis the negative of the outer normal of the
fluid domain Ω,nw=−n. Let the vectors ~eland ~edindicate the directions
of lift and drag. Then, we are able to calculate the lift and the drag with
the boundary integral (1.0.1) where ~elor ~edhave to be inserted instead of ~e.
Here, ~edis the normalized vector directed opposite to the gravity, and ~elis
the normalized vector in the direction opposing the main flow field.
The optimization problem is then formulated as follows: Find a control
uin L2(Γ)nthat maximizes the lift
−ZΓwν∂u
∂n−pn·~elds(1.0.2)
subject to the steady state Navier-Stokes equations describing the motion of
the fluid −ν∆u+ (u·∇)u+∇p= 0 in Ω
div u= 0 in Ω
u=gon Γc,
u= 0 on Γ\Γc,
(1.0.3)
the convex control constraints
g(x)∈Ga.e. on Γc,(1.0.4)
and the maximal drag constraint
−ZΓwν∂u
∂n−pn·~edds≤D0.(1.0.5)
The boundary Γwis a curve satisfying
ZΓw
nds= 0.(1.0.6)
As shown later, the pressure is only unique up to a constant. The constraint
(1.0.6) avoids that this constant changes the objective functional arbitrar-
ily. The control acts on a part of the boundary of the body Γc⊂Γwand
homogeneous Dirichlet boundary conditions are prescribed on the boundary
Γ\Γc.
The set of admissible controls Gis a bounded, convex, closed, and non-
empty subset of Rn. Furthermore, we assume 0∈G, which gives us the
option to turn off the control admissible in the optimization problem. For a
more detailed discussion of such convex control constraints, we refer to [110].
3
CHAPTER 1. INTRODUCTION
Let us shortly review available literature on analysis of optimal control
problems for the Navier-Stokes equations. Starting with Abergel and Temam
[2] there is an ever growing list of contributions. Let us only mention the work
by Gunzburger, Hou, and Svobodny [51], Gunzburger and Manservisi [52],
Hinze and Kunisch [54, 55], Kunisch and de los Reyes [28], de los Reyes and
Yousept [30], de los Reyes and Tröltzsch [29], Abergel and Casas [1], Casas
[24, 23], Tröltzsch and Roubiček [84] and Wachsmuth [103]. Finite-element
error estimates can be found in the work of Casas, Mateos and Raymond
[20]. Optimal flow control problems with state constraints were studied by
Griesse and Reyes [50], Reyes and Kunisch [80].
The novelty of the first part of this thesis is that it combines the use of very
low regular boundary controls, i.e. in L2(Γ), and integral state constraints.
There are only a few contributions to optimal control theory using Dirichlet
controls in L2, see for instance Kunisch and Vexler [61] and Casas, Mateos and
Raymond [20]. In the context of steady-state Navier-Stokes equations this
is a new and promising approach, since the use of L2-controls yields localiz-
able optimality conditions, whereas the use of, for instance, H1/2(Γ)-controls
yields optimality conditions containing non-local boundary operators.
In view of the low L2-regularity of the controls, the boundary integrals
(1.0.2) and (1.0.5) are no longer well-defined, since the velocity field uis not
regular enough to admit traces on the boundary. Therefore, we transform the
boundary integrals into volume integrals leading in case of the drag constraint
to a non-standard mixed control-state constraint, see Section 3.2 below.
As it is well-known, the steady-state Navier-Stokes equations are solvable
in suitable spaces. If the data and/or the Reynolds number 1/ν are small
enough then the solution will be unique. To judge whether this condition is
fulfilled in a concrete application is a delicate issue in particular in the case of
inhomogeneous boundary conditions, see the discussion in the monograph of
Galdi [44]. Hence, instead of assuming smallness of the data, we assume non-
singularity at the optimal control, which is equivalent to unique solvability
of a certain linearized equation, see Section 2.4.2.
By assuming the existence of a linearized Slater point to the state con-
straint (1.0.5), we are able to prove first-order necessary optimality condi-
tions, see Section 4.1. For the special case of smooth controls, the resulting
optimality system simplifies considerably, see Section 4.3. Furthermore, we
state a second-order sufficient optimality condition for the problem under
consideration. The first part of this thesis is complemented by numerical
experiments on a high-lift configuration. One numerical approach was to
solve the associated nonsmooth first-order optimality system of two coupled
Navier-Stokes equations. At the end of this topic, we also implemented an
SQP method with a penalty term in the cost functional to handle the integral
4
state constraint and proved its convergence.
Afterwards, the second part of the paper deals with a nonstationary prob-
lem, considering a model closer to a real setting, that accounts also for turbu-
lence. Here, the flow is computed on the basis of a k−ωWilcox98 turbu-
lence model, see (7.0.2), including the nonstationary Navier-Stokes equations
and high Reynolds numbers. There are many approaches of turbulence mod-
eling. Let us just mention Reynolds-averaged Navier-Stokes (RANS) based
models, one equation models and also two equation models, for instance the
k−model and the k−ωmodels. For every group, one can find many
variations, see for instance [114].
The curse of dimension and the inherent nonlinearity leads to very large
computing times so that a mathematical optimization analogous to the sta-
tionary case is fairly unrealistic. In [19], a generic high-lift configuration was
investigated and one forward solution of the turbulence model took about
48 hours. In the case of the SCCH configuration, which we consider here,
the computation time was nearly twice that number. We think that model
reduction is a method of choice to avoid such extremely long computation
times. To this aim, many authors have considered proper orthogonal decom-
position (POD), see e.g. [4, 62, 63, 111]. The problem is that, in this case,
we have to insert the POD basis as a Galerkin basis in the full turbulence
model, which is a time consuming task. Thus, we decided to follow an alter-
native approach by Luchtenburg, Noack et al. [66, 74] of building a low-order
dynamical system based on uniform oscillation by parameter identification.
The original full optimization problem consists of the cost functional,
calculating the lift, the nonstationary Navier-Stokes equations, the associated
k−ωWilcox98 turbulence model and of constraints on the control. In the
nonstationary part, we discard for simplicity constraints on the drag. The
control function is considered as a periodic function t7→ g,
g=Bcos(ωt),
where Bis the actuation amplitude and ωis the actuation frequency. In
our problem, we consider only the actuation amplitude Bas the optimiza-
tion parameter. For further work, the actuation amplitude is one possible
additional optimization parameter.
In Chapter 7, we consider our optimization problem in detail with the
turbulence model and introduce the standard POD method.
Unfortunately, the POD method does not identify clear frequencies and
amplitudes. Thus, in Chapter 8 we filter the POD mode coefficients in a
way that we neglect fluctuations of frequencies and amplitudes to get clear
5
CHAPTER 1. INTRODUCTION
structures and to calculate the associating modes. With the help of these
coefficients, we build a reduced-order model (ROM), similar to Luchtenburg
et al. [66], consisting only of four ODEs describing uniform oscillations.
We identify the inherent parameters. Based on the results of the ROM, a
formula for the lift calculation is created and so we are able to establish
an optimization problem based on this new lift formula and the low-order
dynamical system.
In Chapter 9, we establish an optimality system, consisting of an objec-
tive cost functional, calculating the lift of the aircraft based on the POD
coefficients, the dynamical system of 4 ODEs as state equations and con-
straints on the control, to investigate the optimization problem numerically
. We need only a fraction of time for one solution of the state equations
compared to the full system while the optimal actuation amplitude is almost
the same.
The main results of the stationary part, except the SQP-method and its
convergence, have been published partially word for word in a joint work with
D. Wachsmuth [58], while Chapter 8 is published partially word for word in
a joint work with B.R. Noack, M. Schlegel, F. Tröltzsch and D. Wachsmuth
[57].
6
Chapter 2
The steady-state Navier-Stokes
equation
The Navier-Stokes equations are a mathematical model to describe the mo-
tion of fluid flow. Claude Louis Marie Henri Navier (b10. February 1785 in
Dijon; d21. August 1836 in Paris) was a French mathematician and physicist
and Sir George Gabriel Stokes (b13. August 1819 in Skreen, County Sligo;
d1. February 1903 in Cambridge) was an Irish mathematician and physicist.
They were the first trying to derive equations of motions for fluid.
The nonlinear incompressible steady state Navier-Stokes system with in-
homogeneous Dirichlet boundary condition is given in its dimensionless form
as follows −ν∆u+ (u·∇)u+∇p=fin Ω
div u= 0 in Ω
u=gon Γc,
u= 0 on Γ\Γc.
(2.0.1)
Here, the velocity field is denoted by uand the pressure by p. In this section,
we investigate the steady-state Navier-Stokes equations mathematically, this
includes the spaces for the solutions. After that, we have to clarify in what
sense the solutions are defined.
2.1 Function spaces
Let us first define the spaces of p-integrable functions and summarize some of
their basic properties. The functions are defined on a domain Ω⊂Rnwith
boundary Γ, which is assumed to be sufficiently smooth, I’ll explain later.
Definition 2.1. Let Ωbe an open subset of Rnand 1≤p < ∞. The set of
7
CHAPTER 2. THE STEADY-STATE NAVIER-STOKES EQUATION
p-integrable functions is defined as
Lp(Ω) = {u: Ω →R;uis measurable and Z
Ω
|u|pdx < ∞},
endowed with the norm
kukLp={Z
Ω
|u(x)|pdx}1/p.
For p=∞, we define
L∞(Ω) = {u: Ω →R;uis measurable and |u(x)| ≤ Ca.e. in Ω,C > 0}
and introduce the norm
kukL∞= inf{C:|u(x)| ≤ Ca.e. in Ω}.
We will provide some basic inequalities to deal with Lebesgue integrable
functions. At first, we have the well-known theorem:
Theorem 2.2 (Hölder).Let u∈Lp(Ω) and v∈Lq(Ω) with 1< p, q < ∞
and 1
p+1
q= 1.
Then uv ∈L1(Ω) and
kuvkL1(Ω) ≤ kukLp(Ω)kvkLq(Ω).
The spaces Lp(Ω) are Banach spaces for 1≤p≤ ∞ and reflexive for
1<p<∞. In L2(Ω),a scalar product can be defined by
(u, v)L2(Ω) =Z
Ω
uv dx
and a Hilbert space structure is also obtained. The space of infinitely dif-
ferentiable functions with compact support is denoted by D(Ω) and its dual,
the distributions space, by D0(Ω).
The Sobolev space Wm,p(Ω) is the space of Lp(Ω) functions whose weak
derivative up to order mis also in Lp(Ω): For these spaces a norm is intro-
duced in the following way:
|u|Wm,p :=
X
|j|<m kDjukp
Lp
1/p
.
8
2.1. FUNCTION SPACES
In the case p= 2, the space Hm(Ω) := Wm,2(Ω) is a Hilbert space with scalar
product
(u, v)Hm=X
|j|≤m
(Dju, Djv)L2.
The closure of D(Ω) in the Wm,p(Ω) norm is denoted by Wm,p
0(Ω). For more
details of Sobolev spaces and a proof for the two following Sobolev imbedding
Theorems, we refer to [3].
Theorem 2.3 (Sobolev imbedding Theorem).Let Ωbe a domain in Rnwith
the cone property. Then the imbeddings
Wm,p(Ω) ,→Lq(Ω),for mp < n and 1≤q < np
n−mp,
Wm,p(Ω) ,→Lq(Ω),for mp =nand 1≤q < ∞,
Wm,p(Ω) ,→C(¯
Ω) ,for mp > n
are continuous.
In the two-dimensional case the imbeddings H1(Ω) ,→Lq(Ω) for 1≤q <
∞and W1,p(Ω) ,→C(¯
Ω) for p > 2are continuous.
Theorem 2.4. Let Ωbe a domain in Rnwith C1-boundary.
1. Suppose mp < n and n−mp < n. This leads to
Wj+m,p(Ω) ,→Wm,q(Ω)
for p≤q≤np/(n−mp).
2. Suppose mp =n. Then, we obtain
Wj+m,p(Ω) ,→Wj,q(Ω)
for p≤q < ∞.
Inhomogeneous boundary values will be defined by the trace of Wk,p(Ω)
on the boundary.
Theorem 2.5. Let Ωbe a bounded Lipschitz-domain and 1≤p≤ ∞. Then
there exists a linear and continuous mapping τ:W1,p(Ω) →Lp(Γ) with
(τu)(x) = u(x)a.e. on Γ
for all u∈W1,p(Ω) ∩C(¯
Ω).
9
CHAPTER 2. THE STEADY-STATE NAVIER-STOKES EQUATION
For a proof of this theorem and the next two, we refer to [3].
Definition 2.6. The element τu is defined by the trace of uon Γand the
mapping is called the trace operator.
Theorem 2.7. Let m≥1,m∈Zand Γof class Cm−1,1. Then, we obtain
that the trace operator τis for mp < n continuous from Wm,p(Ω) to Lr(Γ),
if 1≤r≤(n−1)p
n−mp . For mp =n, we get that τis continuous for all
1≤r < ∞.
Another statement of [3] is:
Theorem 2.8. Let Ωbe of class Cm,m≥1,k∈Z, and 1< p < ∞. Then
the trace operator τis continuous from Wm,p(Ω) to Wm−1/p,p(Γ).
Because we have to deal with functions satisfying div u= 0, we introduce
the space
V:= {u∈ D(Ω)n: div u= 0}
The closure of Vin the H1
0-norm is denoted by Vand if Ωis an open bounded
Lipschitz set, it can be characterized as
V={u∈H1
0(Ω) : div u= 0}.
We assume that the boundary of Ωis in C2. The outer unit normal on Γ
is denoted by n. The boundary Γis the union of mconnected components,
Γ = Sm
j=1 Γj.
Furthermore, we define the following spaces on Ωand Γ:
Hs(Ω) = {v∈Hs(Ω)n: div v= 0 on Ω,hu·n,1iH−1/2(Γj),H1/2(Γj)= 0
∀j∈ {1, . . . , m}}, s ≥0,
Hs
0(Ω) = {v∈Hs(Ω)n: div v= 0 on Ω, u = 0 on Γ}, s ≥1/2,
Hs(Γ) = {v∈Hs(Γ) : RΓju·n= 0 ∀j= 1 . . . m}, s ≥0,
Lp(Ω) = {v∈Lp(Ω)n: div v= 0}, p ≥1,
Lp(Γ) = Lp(Γ)n, p ≥1,
H−s(Γ) = (Hs(Γ))0, s ≥0,
H−s(Ω) = (Hs(Ω) ∩H1
0(Ω))0, s ≥1
Wm,p(Ω) = Wm,p(Ω)n,
Wm,p(Γ) = Wm,p(Γ)n,
10
2.2. THE STOKES EQUATIONS
and the differential operators for vector-valued functions uand scalar-valued
functions p:
∆u∈Rn: (∆u)i= ∆ui, i = 1, . . . , n,
∇p∈Rn: (∇p)i=∂
∂xi
p, i = 1, . . . , n,
∂
∂nu=n∇u∈Rn: ( ∂
∂nu)i=∂
∂nui, i = 1, . . . , n,
(u·∇)u∈Rn: ((u·∇)u)i=
n
X
j=1
uj
∂ui
∂xj
, i = 1, . . . , n.
Remark 2.9. As mentioned later, the terms hu·n,1iH−1/2(Γj),H1/2(Γj)=
0∀j∈ {1, . . . , m}are important to obtain the existence of very weak solutions
for arbitrary large data. The trace operator is only defined for s≥1/2.
Before considering the solvability of the Navier-Stokes equations, we want
to have a look on the Stokes equations.
2.2 The Stokes equations
The Stokes equations are similar to the Navier-Stokes ones, but without the
nonlinear term. They are given by
−ν∆u+∇p=fin Ω
div u= 0 in Ω
u=gon Γc,
u= 0 on Γ\Γc.
(2.2.1)
Let us define the bilinear form a(u, v) := ν(∇u, ∇v). Then we call ua weak
solution of (2.2.1), if
a(u, v) = (f, v)∀v∈V
is satisfied. The next theorem guarantees existence and uniqueness of solu-
tions for (2.2.1) and we refer to [27, 47, 97] for the theory.
Theorem 2.10. For every f∈H−1(Ω) and g∈H1/2(Γ) with
ZΓc
g·nds= 0
the Stokes equation (2.2.1) has a unique solution (u, p)∈H1(Ω) ×L2(Ω),
where pis unique up to a constant.
11
CHAPTER 2. THE STEADY-STATE NAVIER-STOKES EQUATION
The next theorem,see [44, Vol. 1, Thm. IV.6.1], shows that the solution
(u, p)is more regular if the data (f, g)is more regular.
Theorem 2.11. Let ube a solution of the Stokes problem (2.2.1) on a
bounded domain Ω⊂Rn,n≥2of class Cm+2,m≥0corresponding to
f∈Wm,q(Ω), g ∈Wm+2−1/q,q(Γ).
Then, we obtain
u∈Wm+2,q(Ω), p ∈Wm+1,q(Ω).
Moreover, the following estimate holds
kukWm+2,q(Ω) +kpkWm+1,q(Ω) ≤c(kfkWm,q(Ω) +kgkWm+2−1/q,q(Γ))(2.2.2)
with a constant c=c(m, n, q, Ω).
Corollary 2.12. The special case q= 2 and m= 0: For f∈L2(Ω) and
g∈H3/2(Γ), we obtain u∈H2(Ω) and p∈H1(Ω) as solutions for the Stokes
equations (2.2.1). Let us define the associated control-to-state operators
G:H3/2(Γ) →H2(Ω), g 7→ u
with f= 0 and
S:L2(Ω) →H2(Ω), f 7→ u
with g= 0.
We even get the following theorem, see [26, Vol. 6, Thm. 1.11].
Theorem 2.13. Let Ωbe an open bounded set of Rn,n= 2,3, of class Cr,
r= max{m+ 2,2},m≥ −1,m∈Z. Let
f∈Wm,q(Ω) g∈Wm+1,q(Γ)
satisfy the compatibility condition
Z
Γc
g·nds= 0.
Then there exists a unique solution
(u, p)∈Wm+2,q(Ω) ×Wm+1,q(Ω)
of the Stokes problem (2.2.1) (pis unique up to a additive constant), satisfying
(2.2.2).
Let us specify, in what sense we want to solve the Navier-Stokes equations.
The first and standard way is to consider weak solutions.
12
2.3. WEAK FORMULATION
2.3 Weak formulation
See [44] for more details on weak solutions of (1.0.3). For simplicity, we
define some functions b:H1(Ω) ×H1(Ω) ×H1(Ω) →R
b(u, v, w) := ((u·∇)v, w)L2(Ω) (2.3.1)
and B:H1(Ω) 7→ (H1(Ω))0for u, v ∈H1(Ω)
hB(u), vi(H1(Ω))0,H1(Ω) := b(u, u, v).(2.3.2)
Due to to the quadratic character of B, its differentiability is easily to see.
The first Fréchet derivative B0(¯u)uof Bwith respect to ¯uis a functional of
(H1(Ω))0and has the following form
hB0(¯u)u, vi(H1(Ω))0,H1(Ω) =b(¯u, u, v) + b(u, ¯u, v)
with v∈H1(Ω). The second Fréchet derivative B00(¯u)[˜u, ˆu]is given by
hB00(¯u)[˜u, ˆu], vi(H1(Ω))0,H1(Ω) =b(˜u, ˆu, v) + b(ˆu, ˜u, v).
Additionally, we get the following properties for the nonlinearity band B,
see [55],
Lemma 2.14. The nonlinear term band Bfulfills
1. |b(u, v, w)| ≤ ckuk1/2
L2(Ω)kuk1/2
H2(Ω)kvkH1(Ω)kwkL2(Ω),∀u∈H2(Ω), v ∈
H1(Ω), w ∈L2(Ω),
2. |b(u, v, w)| ≤ ckukH1(Ω)kvkH1(Ω)kwkH1(Ω),∀u, v, w ∈H1(Ω),
3. hB00(¯u)[˜u, ˆu], vi(H1(Ω))0,H1(Ω) =hB0(˜u)ˆu, vi(H1(Ω))0,H1(Ω)
4. 1
2hB00(¯u)[u, u], vi(H1(Ω))0,H1(Ω) =hB(u), vi(H1(Ω))0,H1(Ω) for u= ˜u= ˆu
with a constant c∈R.
Multiplying (1.0.3) by test functions (v, q)∈V×L2(Ω), we obtain by
partial integration the following weak formulation
a(u, v) + b(u, u, v) = hf, viV0,V in Ω(2.3.3)
τu =gon Γc(2.3.4)
where τ:H1:→H1
2(Γc)is the trace operator.
13
CHAPTER 2. THE STEADY-STATE NAVIER-STOKES EQUATION
Theorem 2.15. Let Ω∈Rnbe a bounded locally Lipschitz domain of Rn,n=
2,3with Γ := ∂Ωconstituted by m+ 1 connected components Γ1,...,Γm+1,
m≥0,
f∈H−1(Ω) and g∈H1/2(Γc)
with ZΓc
g·nds= 0.
Then, there exists under additional assumptions, see [44, Vol. 2, Thm. VIII.4.1,
Thm. VIII.4.2], at least one solution (u, p)∈H1(Ω) ×L2(Ω) for (2.0.1) sat-
isfying the estimate
kpkL2(Ω) ≤c(kfkH−1(Ω) +kuk2
H1(Ω) +νkukH1(Ω)).
Furthermore, there exists c1=c1(n, Ω) such that if
kgk1/2,2(Γ) ≤c1ν/2,
uverifies
kukH1(Ω) ≤c2(kfkH−1(Ω) +kgk2
H1
2(Γ) + (1 + ν)kgkH1
2(Γ)).
with c2=c2(n, Ω). If
kfkH−1(Ω) +kgk2
H1
2(Γ) + (1 + ν)kgkH1
2(Γ) < c3ν2
is additionally satisfied with c3= min{c1,1/c2k}, then uis unique and pis
unique up to a constant.
Proof. For the proof and the additional assumptions, we refer to [44, Vol. 2,
Thm. VIII.4.1, Thm. VIII.4.2].
We can obtain some extra regularity of the solution (u, p), if the right
hand side is smooth enough, see [44, Vol. 2, Thm. VIII.5.2]:
Theorem 2.16. Let Ωbe a bounded domain of Rn,n≥2, of class C2. Let
u∈W1,2(Ω) ∩Ln(Ω)
be divergence-free, satisfying (2.3.3) for some p∈L2(Ω) and for all v∈
C∞
0(Ω). Then, if
f∈Lq(Ω), g ∈W2−1/q,q(Γc)
with
q∈(1,∞),if n= 2,
14
2.4. VERY WEAK FORMULATION
while
q∈(2n/(n+ 2),∞),if n > 2,
it follows that
(u, p)∈W2,q(Ω) ×W1,q(Ω).
Moreover, if Ωis of class Cm+2 and
f∈Wm,q(Ω), g ∈Wm+2−1/q,q(Γc)
with m≥1and
q∈(1,∞),if n= 2,
while
q∈(n/2,∞),if n > 2,
then
(u, p)∈Wm+2,q(Ω) ×Wm+1,q(Ω).
We obtain by Theorem 2.16 with q= 2:
Corollary 2.17. Let f∈L2(Ω) and g∈H3/2(Γc). Then, it follows that
(u, p)∈H2(Ω) ×H1(Ω).
2.4 Very weak formulation
Here, we have to resort to the notation of very weak solutions, since in general
weak solutions in u∈H1(Ω) do not exist due to the desired low regularity
L2(Γc)of boundary data. The theory in this section is based on [69, 40, 39].
Definition 2.18. Let g∈H0(Γ) be given. Then we call u∈L2n/(n−1)(Ω)
avery weak solution of the state equation (2.0.1) if for all test functions
v∈H2
0(Ω),π∈H1(Ω) it holds
ZΩ
(u·(−ν∆v)−(u·∇)vu)dx+ZΓ
g·ν∂v
∂nds=hf, viH−1(Ω),H1(Ω) (2.4.1a)
and ZΩ∇π·udx−ZΓ
(g·n)πds= 0.(2.4.1b)
15
CHAPTER 2. THE STEADY-STATE NAVIER-STOKES EQUATION
Here, the first equation is obtained by partially integrating the Navier-
Stokes equations twice and using the equation:
ZΩ
(u·∇)u·vdx=
n
X
i,j=1 ZΩ
ui
∂uj
∂xi
vjdx
=
n
X
i,j=1 −ZΩ
ui
∂vj
∂xi
uj+∂ui
∂xi
vjujdx+ZΓ
uiujnivjds
=−ZΩ
(u·∇)v·u+u·vdiv udx+ZΓ
(u·n)u·vds
Then, we can reformulate the first equation of (2.0.1) to (2.4.1a) by
ZΩ
(−ν∆u+ (u·∇)u+∇p)vdx−hf, viH−1(Ω),H1(Ω)
=ZΩ
ν∇u·∇v−(u·∇)v·u+u·vdiv u+p·div vdx
+ZΓ
(u·n)u·v−ν∂u
∂nvds−hf, viH−1(Ω),H1(Ω)
=ZΩ
u·(−ν∆v)−(u·∇)v·udx+ZΓ
g·ν∂v
∂nds−hf, viH−1(Ω),H1(Ω).
The second equation (2.4.1b) is the weak formulation of div u= 0:
ZΩ
div uπdx=−ZΩ
u·∇πdx+ZΓ
(u·n)πds
=−ZΩ
u·∇πdx+ZΓc
(g·n)πds.
Moreover as discussed in [40], it incorporates the Dirichlet boundary condi-
tion for the normal component of u, since the term RΓg·ν∂v
∂ndsacts only
on tangential components. For example, let n= 3 and and vbe a smooth
function. Then it holds ∂v
∂n(curlv)×n, which implies that ∂v
∂nis orthogonal to
the outer normal, and thus the product RΓu·ν∂v
∂ndsin (2.4.1a) acts only on
the tangential component of u.
The existence of very weak solutions with inhomogeneous Dirichlet bound-
ary conditions is discussed in [39, 40, 45, 69]. We remark that it is essential
to have RΓig·nds= 0 for all connected components of Γto obtain existence
of solutions for arbitrary large data. For boundary data in H0(Γ) it holds
the following.
Remark 2.19. If the data are regular and the problem has a variational
solution (u, p)∈H1(Ω) ×L2(Ω), then it is easy to see that the variational
solution is also a very weak solution.
16
2.4. VERY WEAK FORMULATION
In the next theorem, we investigate the linearized state equation, see [69,
Theorem 3].
Theorem 2.20. Let g∈L2(Γc)and z∈L2n/(n−1)(Ω)2. Then the linearized
problem
−ν∆u+ (z·∇)u+∇p=fin Ω
div u= 0 in Ω
u=gon Γc
has a unique very weak solution u∈L2n/(n−1)(Ω) and additionally, we obtain
kukL2n/(n−1)(Ω) < c(1 + kzkL2n/(n−1)(Ω))kgkL2(Γc).
Let us now formulate the main result of this section.
Theorem 2.21. For every f∈H−1(Ω) and g∈H0(Γ), there exists a very
weak solution u∈L2n/(n−1)(Ω) of (2.0.1). In the two-dimensional case, this
solution belongs to H1/2(Ω).
The existence proof and a discussion of the smallness assumption of the
data and/or the Reynolds number Re = 1/ν can be found in [69]. The
H1/2-regularity for the 2d-case result can be proven following the lines of
[69]. Unique solvability with respect to less regular data, i.e. in W−1/q,q(Γ)
is investigated in the articles by Farwig, Galdi, Sohr [39, 40, 45]. Once
existence of a solution is proven, the pressure field pcan be reconstructed by
means of De Rham’s Lemma, see for instance [97].
One can find in [69] that there exists a distribution p∈W−1,2n/(n−1) such
that
−∆u−(u·∇)u+∇p= 0
holds in the sense of distributions.
In view of the existence result, let us define for abbreviation the state
space
U:= L2n/(n−1)(Ω).
2.4.1 More regular solutions
Let us briefly show that more regular boundary data in Lp(Γ) yields more
regular solutions. In the following considerations, we will split the state
uin two parts, u=u0+u1. The function u0is defined as the unique very
weak solution to the Stokes equation with inhomogeneous Dirichlet boundary
conditions −ν∆u0+∇p0= 0 in Ω
div u0= 0 in Ω
u0=gon Γ.
(2.4.2)
17
CHAPTER 2. THE STEADY-STATE NAVIER-STOKES EQUATION
Then u1=u−u0solves the following equation subject to homogeneous
Dirichlet boundary conditions
−ν∆u1+ (u·∇)u1+∇p1=−(u·∇)u0in Ω
div u1= 0 in Ω
u1= 0 on Γ.
(2.4.3)
As one can easily see, both systems are uniquely solvable. At first, let us
prove higher regularity of u0.
Lemma 2.22. Let g∈Lp(Γ) ∩H0(Γ),p≥2, be given. Then the solution of
(2.4.2) satisfies u0∈Lq(Ω), where qis given by
q=(np
n−1if 2≤p < ∞,
+∞if p=∞, n = 3.(2.4.4)
Proof. The mapping g7→ u0is linear and continuous from W−1/p,p(Γ) to
Lp(Ω) and from W1−1/p,p(Γ) to W1,p(Ω), see e.g. [12, 25]. By interpola-
tion arguments, we have continuity of this solution mapping from Lp(Γ) to
W1/p,p(Ω). The claim follows by the imbedding argument W1/p,p(Ω) ,→
Lq(Ω), see Theorem 2.3. The result for p=∞,n= 3 can be found in
[91].
Applying this result, we can prove higher regularity of the function u1
and in consequence of the solution uof the nonlinear system.
Lemma 2.23. If the boundary data gis in Lp(Γ) ∩H0(Γ),p≥2 (n= 2) or
p≥4 (n= 3), then ubelongs to Lq(Ω) with qgiven by (2.4.4).
Proof. Let us first consider the 2d-case, n= 2. Then we have by Theo-
rem 2.21 u∈L4(Ω). This implies (u· ∇)u∈H−1(Ω), hence u1=u−u0
solves −ν∆u1+∇p1=−(u·∇)uin Ω
div u1= 0 in Ω
u1= 0 on Γ,
(2.4.5)
and we have u1∈H1
0(Ω) ,→Lt(Ω) for all 0< t < ∞.
In the 3d-case, n= 3, Theorem 2.21 gives u∈L3(Ω). Since p≥4by
assumption, Lemma 2.22 yields u0∈L3p/2(Ω) ⊂L6(Ω). Taking v∈H1
0(Ω),
we find after integration by parts
|ZΩ
(u·∇)u0vdx|=|−ZΩ
(u·∇)vu0dx| ≤ ck∇vkL2(Ω)ku0kL6(Ω)kukL3(Ω).
18
2.4. VERY WEAK FORMULATION
Then (u·∇)u0belongs to H−1(Ω), and the solution u1=u−u0of (2.4.3)
belongs to H1
0(Ω) ,→L6(Ω). Hence u=u0+u1is in L6(Ω) as well. This in
turn gives (u·∇)u∈W−1,r(Ω),r=6p
p+4 ≥4with the following argument:
|ZΩ
(u·∇)u0vdx| ≤ ck∇vkLr0(Ω)ku0kL3p/2(Ω)kukL6(Ω).
By [25], the function u1as solution of (2.4.5) belongs then to W1,4(Ω) ,→
L∞(Ω).
In both cases, (n= 2,3), the function u1is as regular as u0, hence the
solution ubelongs to the space Lq(Ω) with qas in (2.4.4).
2.4.2 Regularity assumption
It is well known that the stationary Navier-Stokes equations are uniquely
solvable if the data/the control function gis small, see for example [69,
Theorem 4].
Hence, if we want to have a unique response uto each control gwe would
have to impose restrictions on the control to enforce uniqueness of solutions.
This technique is widely employed in optimal control of the stationary Navier-
Stokes equations, see e.g. [50, 80, 81, 84, 103]. We will however proceed
without a smallness assumption and therefore with non-uniqueness of the
solutions. Since we allow multiple solutions of the state equation, we have
to clarify the meaning of optimality.
Definition 2.24. A pair (¯u, ¯g)is called locally optimal, if there exist ρu, ρg>
0such that J(¯u, ¯g)≤J(u, g)holds for all admissible pairs (u, g)with ku−
¯ukU< ρuand kg−¯gkL2(Γc)< ρg.
Here, a pair (u, g)is admissible if it satisfies the constraints (1.0.3)–
(1.0.5).
Instead of enforcing uniqueness of solutions for all controls, we will impose
the following regularity condition on an optimal state. A similar assumption
is used to derive error estimates for distributed control in [20].
Definition 2.25. A pair (¯u, ¯g)∈U×H0(Γ) is called non-singular, if the
linearized Navier-Stokes equation
−ν∆u+B0(¯u)u+∇p=fin Ω
div u= 0 in Ω
u=gon Γ,
(2.4.6)
19
CHAPTER 2. THE STEADY-STATE NAVIER-STOKES EQUATION
admits a unique very weak solution u∈Ufor all g∈H0(Γ) and f∈H−1(Ω).
Moreover, we assume that the solution mapping g7→ ufor f= 0 is linear
and continuous from H0(Γ) to H1/2(Ω), and the mapping f7→ ufor g= 0
is linear and continuous from H−1(Ω) to H1
0(Ω).
This condition is fulfilled, if the state ¯uis small enough [79, Lemma
B.1]. The assumption of non-singularity implies that the state equation can
be solved uniquely in a neighborhood of the reference/optimal control and
state, confer [20, Theorem 2.5] for a proof using an implicit function theorem.
Theorem 2.26. Let (¯u, ¯g)∈U×H0(Γ) be a non-singular solution of (1.0.3).
Then there exist an open neighborhood O(¯g)of ¯gin H0(Γ), an open neigh-
borhood O(¯u)of ¯uin U, and a mapping Sfrom O(¯g)⊂H0(Γ) to O(¯u)⊂
U=L2n/(n−1)(Ω) of class C2such that, for all g∈ O(¯g),S(g) = uis the
unique very weak solution in O(¯u)of (1.0.4).
Furthermore, the action of the first Fréchet derivative u=S0(¯g)gis given
by the unique very weak solution of the linearized equation (2.4.6).
20
Chapter 3
The optimal control problem
Since we will work with very weak solutions u∈U=L2n/(n−1)(Ω), we have
to clarify the meaning of the boundary integrals in the objective functional
(1.0.2) and the state constraint (1.0.3). These would be well-defined if the
regularity ∂u
∂n∈L1(Γ) could be guaranteed. This is not fulfilled for very weak
solutions from H1/2(Ω) or U. We will thus extend the linear functionals
in (1.0.2) and (1.0.3) to the larger space U. The idea of reformulating the
functionals in this way arises by reading the paper of Braack and Richter
[18].
In the first section of this chapter, we introduce reformulated boundary
integrals for (1.0.2) and (1.0.3) with once partial integration. They are well-
defined for u∈H1(Ω) and we need them in Section 4.3. Here, we consider
g∈Gand the associated state u∈U. Therefore, we have to reformulate
(1.0.2) and (1.0.3) again with partially integrating twice in the second section.
3.1 Reformulation of the boundary integrals
Let us assume for a while that u, p ∈C1(Ω)nare classical solutions of (1.0.4)
to the control g. Multiplying the state equations (1.0.4) with a function
ϕi∈H1(Ω), we obtain by partial integration
0 = (−ν∆u+ (u·∇)u+∇p, ϕi)
=ZΩ
(ν∇u·∇ϕi+ (u·∇)uϕi)dx−ZΓν∂u
∂n−pnϕids.
In order to represent the functionals in (1.0.2), (1.0.3), let us introduce func-
tions ϕi,i∈ {d, l}, that take suitable chosen boundary values. Let ϕidenote
21
CHAPTER 3. THE OPTIMAL CONTROL PROBLEM
the weak solutions of the Stokes equations
−∆ϕi+∇πi= 0 in Ω
div ϕi= 0 in Ω(3.1.1)
with boundary values
ϕl=(~elon Γw
0on Γ\Γw
, ϕd=(~edon Γw
0on Γ\Γw
.
Using the transformation above, we can write for i∈ {d, l}.
−ZΓwν∂u
∂n−pn~eids=−ZΓν∂u
∂n−pnϕids
=ZΩ
(ν∇u·∇ϕi+ (u·∇)uϕi)dx.
Taking the right-hand side of this expression, we define the functional
ˆ
F~ei(u) = ZΩ
(ν∇u·∇ϕi+ (u·∇)uϕi)dx(3.1.2)
for i∈ {d, l}. This function is well-defined for u∈H1(Ω), but for u∈Uwe
have to reformulate the boundary integrals once again.
3.2 Further reformulation of the boundary in-
tegrals
Let us now assume that u∈C2(Ω)n, p ∈C1(Ω)nare a classical solution of
(1.0.4) to g. Multiplying the state equation (1.0.4) with ϕi∈H2(Ω), we
obtain by partially integrating twice
0 = (−ν∆u+ (u·∇)u+∇p, ϕi)
=ZΩ
(ν∇u·∇ϕi+ (u·∇)uϕi)dx−ZΓν∂u
∂n−pnϕids
=ZΩ
(−νu ·∆ϕi−(u·∇)ϕiu)dx
+ZΓνu∂ϕi
∂n+ (u·n)(ϕi·u)ds−ZΓν∂u
∂n−pnϕids,
22
3.2. FURTHER REFORMULATION OF THE BOUNDARY
INTEGRALS
where ϕi,i∈ {d, l}are defined as in the section before. Here, we need
ϕi∈H2(Ω),i∈ {d, l}. Using the transformation above, we can write for
i∈ {d, l}.
−ZΓwν∂u
∂n−pn·~eids=−ZΓν∂u
∂n−pnϕids
=ZΩ
(νu ·∆ϕi+ (u·∇)ϕiu)dx
−ZΓνu∂ϕi
∂n+ (u·n)(ϕi·u)ds.
Taking the right-hand side of this expression, we define the functional
˜
Fi(u) = ZΩ
(νu ·∆ϕi+ (u·∇)ϕiu)dx−ZΓνu∂ϕi
∂n+ (u·n)(ϕi·u)ds.
This function is well-defined for u∈Hs(Ω),s > 1/2, but not for u∈U. To
handle this problem, we substitute the state uby the control function u=g
on the boundary. This yields
˜
F~ei(u, g) = Z
Ω
(νu ·∆ϕi+ (u·∇)ϕiu)dx−Z
Γνg∂ϕi
∂n+ (g·n)(ϕi·g)ds.
(3.2.1)
In contrast to (1.0.2) and (1.0.3), the function ˜
F~ei,i∈ {d, l}, is well-defined
for states u∈Uand controls g∈H0(Γ), since the functions ϕare very
regular in comparison to the very weak solutions in U. Their boundary values
are in fact a constant vector, thus, the regularity of ϕiis only influenced by
the regularity of the boundary Γ. But we have to note that the control now
appears nonlinearly in the functionals which leads to additional difficulties
for the optimal control problem, see Subsection 3.3.1.
The following result can be deduced from Theorem 2.11:
Lemma 3.1. The functions ϕi,i∈ {d, l}, belong to H2(Ω) ∩W2,p(Ω) for all
p < ∞.
Then we can prove a continuity and differentiability statement for fi.
Lemma 3.2. The functions ˜
F~eiare continuous and twice Fréchet- differen-
tiable from U×H0(Γ) to R. Moreover, for given u∈Lq(Ω) and g∈Lp(Γ),
p, q ∈(1,+∞), it holds ˜
F0
ei,u(u, g)∈Lq(Ω) and ˜
F0
ei,g(u, g)∈Lp(Ω).
Proof. We get the continuity and the Fréchet-derivatives by standard argu-
ments. Because of Lemma 3.1, ϕi,i∈ {d, l}possesses enough regularity to
give ˜
F0
ei,u(u, g)∈Lq(Ω) and ˜
F0
ei,g(u, g)∈Lp(Ω).
23
CHAPTER 3. THE OPTIMAL CONTROL PROBLEM
3.3 The optimal control problem
Now, we are able to reformulate the original optimal control problem (1.0.2)–
(1.0.5) using the very weak solutions and the extended functionals ˜
F~el,˜
F~ed
(3.2.1). Let us denote the following optimal control problem (3.3.1)–(3.3.4)
by (P): Minimize
J(u, g) := −˜
F~el(u, g) + α
2kgk2
L2(Γc)(3.3.1)
subject to the very weak form of
−ν∆u+ (u·∇)u+∇p=fin Ω
div u= 0 in Ω
u=gon Γc,
u= 0 on Γ\Γc,
(3.3.2)
the control constraints
g∈Gad := {v∈H0(Γc) : (ga)i(x)≤vi(x)≤(gb)i(x)
a.e. on Γc,∀i∈ {1, . . . , n}},(3.3.3)
and the integral control-state constraint
˜
F~ed(u, g)≤D0.(3.3.4)
The functions gaand gbare elements of L∞(Γc)with (ga)i≤(gb)ifor
all i∈ {1, . . . , n}a.e. on Γcand we assume 0∈Gad.Here, we introduced
an additional regularization term α
2kgk2
L2(Γc), where αis called the Tichonov
parameter, which measures the cost of the control and is important for the
optimality system, see Chapter 4. The parameter αis supposed to be posi-
tive.
3.3.1 Existence of solutions
After introducing the optimality problem, we would like to prove the exis-
tence of solutions of problem (P). Unfortunately, we can not prove that the
objective functional is bounded from below. This is due to the absence of a
uniqueness result for the state equation for large data. Moreover, bounds on
the state of the kind kukU≤CkgkL2(Γc)can only be derived for small data.
This is different to the distributed control problems for Navier-Stokes equa-
tions, where we can test the state equation with the state itself and obtain
an a-priori bound without smallness assumptions.
24
3.3. THE OPTIMAL CONTROL PROBLEM
Due to the fact that we are not able to prove the existence of solutions
for the problem (P), we will introduce a modification. Let us consider the
minimization of
˜
J(u, g) := ψ(−˜
F~el(u, g)) + ˜α
2kuk2
H1/2(Ω) +α
2kgk2
L2(Γc).(3.3.5)
Here, ˜αand αare positive and small parameters. The function ψ:R→Ris
assumed to be continuous, monotone increasing and bounded from below; e.g.
ψ(r)≥ψmin for all r∈R. For example, one can choose the function ψ(r) =
−log(L0−r)with L0∈R. This function additionally forbids situations
where the lift is too small; that is, smaller than a prescribed value L0, because
(L0−r)∈R+has to be fulfilled.
It appears that the functionals ˜
F~eiare not weakly continuous with respect
to g∈H0(Γc). In order to overcome this difficulty, we will impose the
following control constraint, where ˜
Gad is a closed and convex set such that
˜
Gad ⊂g∈Gad :ZΓw
(g·n)(ϕi·g)ds= 0, i ∈ {d, l}.(3.3.6)
If the control boundary is not part of the observation boundary, i.e. Γc∩Γw=
∅, one can choose ˜
Gad =Gad. This choice is also valid in the case of pure
tangential controls g·n= 0.
Let us denote the modified minimization problem (3.3.5)–(3.3.6) by (˜
P).
Theorem 3.3. If there is an admissible pair (u0, g0)∈H1/2(Ω)ט
Gad, which
satisfies (3.3.2)–(3.3.4) and the control constraint (3.3.6), then the problem
(˜
P) admits at least one solution.
Proof. The objective functional ˜
Jis bounded from below by construction.
We can restrict the optimization problem to the set of all admissible pairs
(u, g)with ˜
J(u, g)≤˜
J(u0, g0)without changing the set of global minimizers.
Let us take such an admissible pair. We then obtain
˜α
2kuk2
H1/2(Ω) +α
2kgk2
L2(Γc)≤ −ψ(−˜
F~el(u, g)) + J(u0, g0)
≤ −ψmin +J(u0, g0),
(3.3.7)
which implies that the set of admissible pairs with smaller value of the ob-
jective than J(u0, g0)is bounded.
Since ˜
Jis bounded from below, there exists a minimizing sequence
(un, gn)∈H1/2(Ω)×H0(Γ). In view of (3.3.7), this sequence is bounded and
we can extract a weakly converging subsequence, which is again denoted by
(un, gn), i.e. un*¯uin H1/2(Ω) and gn*¯gin H0(Γ).
25
CHAPTER 3. THE OPTIMAL CONTROL PROBLEM
By compact embeddings, we have un→uin Lp(Ω)nfor all p < 3after
extracting another subsequence, see [3, Theorem 6.2].
The set ˜
Gad is weakly closed by construction, ˜
Gad is defined as closed and
convex, which implies ¯g∈˜
Gad. Together with the control constraint gn∈
˜
Gad, this allows us to pass to the limit in the functions ˜
F~ei,lim ˜
F~ei(un, gn) =
˜
F~ei(¯u, ¯g).
Passing to the limit in the very weak solution is straight-forward, which
implies that ¯uis a very weak solution to ¯g. Now, standard arguments using
lower semi-continuity of norms conclude the proof.
Let us summarize the difficulties in proving the existence of solutions:
1. The functional (3.3.1) is not bounded from below, since there is no
a-priori boundary kukU≤CkgkL2(Γc).
2. If a minimizing sequence would exist, the sequence unis not necessarily
bounded in U.
3. The functions ˜
F~eiare not weakly continuous on U×H0(Γ).
Therefore the modification of the objective and the control constraint were
made to cope with these points.
1. The function ψguarantees that the objective function is bounded from
below.
2. The regularization term kuk2
H1/2(Ω) gives boundedness of unin H1/2(Ω).
By compact embeddings, it allows furthermore to pass to the limit in
the part of ˜
F~eithat involves the state u.
3. The control constraint RΓ(g·n)(ϕi·g)ds= 0 permits to pass to the
limit in the non-linear part of ˜
F~eithat involves the control.
Of course, there are several other possibilities to enforce existence of solu-
tions. For instance, one could add regularization with respect to stronger
norms in g. This would, however, lead to different first-order necessary op-
timality conditions, which are more challenging to solve numerically; see the
comments below.
Another popular change would be to explicitly impose a smallness con-
dition on the controls. However, known smallness conditions are difficult to
verify, especially in the case of non-homogenous Dirichlet boundary condi-
tions when Γconsists of more than one connected component, see [44]. But
this is the case in the application we have in mind. Moreover, the smallness
26
3.3. THE OPTIMAL CONTROL PROBLEM
condition on 1/ν =Re is expected to be violated for concrete applications,
where typically 1/ν =Re is large.
In the following, we consider the original problem without the function
ψand the term ˜α
2kukH1,2(Ω). These modifications, particularly the term
with the norm of H1,2(Ω), would generate problems with the variational
inequality, for example in the Subsection 4.4, where we consider the second-
order sufficient optimality conditions.
27
CHAPTER 3. THE OPTIMAL CONTROL PROBLEM
28
Chapter 4
Optimality conditions
In this chapter, we will establish the first-order necessary and the second-
order sufficient optimality conditions. These conditions are very important
for the numerical solution methods in the next chapter. First, we will inves-
tigate the first-order necessary optimality condition for the solution in the
very weak sense. We will also introduce a new space for the control function,
fitting better to the numerical investigation.
4.1 First order necessary optimality conditions
Let us return to problem (P) stated at the beginning of Section 3.3. We
will now derive necessary optimality conditions to characterize local optimal
solutions. Here, we will follow the presentation in [102, Section 6.1.2] of the
regularity theory of [115]. Let now (¯u, ¯g)∈U×Gad be a non-singular and
locally optimal pair for (P). and we define the operator
G= (G1, G2) : O(¯u)×O(¯g)→U×R
(u, g)7→ S(g)−u
˜
F~ed(u, g)−D0(4.1.1)
and the cone
K={0}×(−∞,0] ⊂U×R.
Here, D0denotes the drag constraint (1.0.5) and Kinduces a partial ordering
<Kon U×Rby: x <K0⇔x∈K. It is easy to show using Theorem 2.26
that the mapping Gis twice Fréchet differentiable. Then, (¯u, ¯g)is a local
solution of the minimization problem
min
u∈O(¯u),g∈O(¯g)J(u, g)subject to G(u, g)<K0, u ∈Gad.
29
CHAPTER 4. OPTIMALITY CONDITIONS
Let us define the Lagrangian associated with this optimization problem
L(u, g;θ, ξ) = J(u, g)+hS(g)−u, θiU0,U+ξ(˜
F~ed(u, g)−D0), θ ∈U0, ξ ∈R.
To show existence of Lagrange multipliers, we will assume the following reg-
ularity condition: there exists ˜g∈Gad ∩ O(¯g)and the associated state
˜u∈ O(¯u)such that
G0
1(¯u, ¯g)(˜u, ˜g) = 0, G2(¯u, ¯g) + G0
2(¯u, ¯g)(˜u−¯u, ˜g−¯g)<0.(4.1.2)
This condition is sufficient for the Zowe-Kurcyusz regularity assumption
[115], see e.g. [99]. Now, we are able to establish the necessary first-order
optimality conditions for (P).
The optimality system
Under these assumptions, the existence of Lagrange multipliers follows by
known results, see e.g. [99, 102, 115].
Theorem 4.1. Let (¯u, ¯g)a non-singular local optimal solution for (P). Let
us assume that there are ˜g∈Gad ∩O(¯g)and the associating state ˜u∈ O(¯u)
such that the linearized Slater condition (4.1.2) is satisfied. Then there exists
θ∈U0and ξ≥0such that the equation
θ=−˜
F0
~el,u(¯u, ¯g) + ξ˜
F0
~ed,u(¯u, ¯g),(4.1.3a)
the variational inequality
(α¯g−˜
F0
~el,g(¯u, ¯g)+ξ˜
F0
~ed,g(¯u, ¯g)+S0(¯g)∗θ, g−¯g)L2(Γc)≥0∀g∈Gad,(4.1.3b)
and the complementarity condition
ξ(˜
F~ed(¯u, ¯g)−D0)=0, ξ ≥0,˜
F~ed(¯u, ¯g)≤D0(4.1.3c)
hold.
Here, the adjoint operator S0(¯g)∗:U0→H0(Γ)0appears with
U0=L(2n/(n−1))0(Ω) = L2n/(n+1)(Ω).
The dual space of Lp(Ω) can be identified with the space Lp0(Ω), where p0is
defined by p0:= p/(p−1). It is related to the solution of the so-called adjoint
equation. In fact, it holds [79]:
30
4.1. FIRST ORDER NECESSARY OPTIMALITY CONDITIONS
Theorem 4.2. The action of S0(¯g)∗can be characterized as follows: for
given θ∈U0it holds
S0(¯g)∗θ=−ν∂λ
∂n−πnΓc
,(4.1.4)
where λ∈H1
0(Ω) ∩W2,r(Ω)nis the unique solution of the equation in a weak
sense
ZΩ
(ν∇λ·∇v−(¯u·∇)λv −(v·∇)λ¯u)dx=hθ, viH−1,H1(4.1.5)
for all v∈H1
0(Ω) and π∈W1,r(Ω) the associated pressure field for all
r∈[2,∞).
Proof. The representation of S0(¯g)∗is proven for instance in [79]. It remains
to investigate the regularity of λ. The right-hand side hθ, viH−1,H1is given
according to the previous Theorem 4.2 by
hθ, viH−1,H1=ZΩ
(−ν∆ ˜ϕv −(v·∇) ˜ϕ¯u−(¯u·∇) ˜ϕv)dx,
where we used the notation ˜ϕ=−ϕl+ξϕd. By assumption, Gad is a subset
of L∞(Γ), hence ¯g∈L∞(Γ) and ¯u∈Lq(Ω),2≤q < ∞for n= 2,2≤q≤ ∞
for n= 3, cf. (2.4.4). Due to the high regularity of ˜ϕ, compare Lemma 3.1,
we can estimate with some ˜p>nsuch that W2,˜p(Ω) ,→W1,∞(Ω)
|hθ, viH−1,H1| ≤ c(1 + k¯ukLq)k˜ϕkW2,˜pkvkLr
with 1/q + 1/˜p+ 1/r = 1. Since qand ˜pcan be chosen arbitrarily large
(but not equal to ∞), we obtain θ∈Lq(Ω)n, for all q∈(2,∞). Let us now
estimate the addend on the left-hand side of (4.1.5) that comes from the
nonlinearity of the state equation:
ZΩ
((¯u·∇)λv + (v·∇)λ¯u)dx≤ck¯ykLqk∇λkLpkvkLr(4.1.6)
with 1/q + 1/p + 1/r = 1,2≤q < ∞.
Since (¯u, ¯g)is non-singular, the equation (4.1.5) is uniquely solvable with
solution λ∈H1
0(Ω). That is, estimate (4.1.6) holds with p= 2. We can
interprete the adjoint state as the weak solution of a Stokes equation, where
the terms in (4.1.6) are put on the right-hand side. This allows to apply
known regularity results for the Stokes equation.
With p= 2 the estimate (4.1.6) holds for all r > 2, hence the functional in
(4.1.6) is in Lr0(Ω) for all r0<2. The regularity result by Galdi [44, Lemma
31
CHAPTER 4. OPTIMALITY CONDITIONS
IV.6.1] gives in a first step the regularities λ∈W2,r0(Ω) and π∈W1,r0(Ω)
for all r0<2.
By embedding arguments, we have then ∇λ∈Lp(Ω)n,n, where pdepends
on n:p∈(2,∞)for n= 2;p∈(2,6) for n= 3.
In the 2d-case, (4.1.6) is valid for all p < ∞, which allows us to chose r
arbitrary small with r > 1. That is, the functional involving ¯yand ∇λis
in Lq(Ω) for all q < ∞. Again applying the regularity result for the Stokes
equation, we find λ∈W2,q(Ω) and π∈W1,q(Ω) for all q < ∞.
For the three-dimensional case, we obtain similarly λ∈W2,q(Ω) and π∈
W1,q(Ω) for all q < 6. By continuous imbeddings, ∇λ∈L∞(Ω)n,n follows,
which gives after applying again Galdi’s regularity result λ∈W2,q(Ω) and
π∈W1,q(Ω) for all q < ∞.
The adjoint pressure πis only determined up to constant. This fact is
usually circumvented by requiring RΩπdx= 0. Here, it is not necessary to
fix the constant, since the constant does not change the variational inequality
due to the construction of H0(Γ)
(πn, v)L2(Γc)= ((π+c)n, v)L2(Γc)∀c∈R, v ∈H0(Γ).
Furthermore, the variational inequality (4.1.3b) can be written as a non-
smooth equation.
Corollary 4.3. Let the assumptions of the previous theorem be fulfilled.
Then the variational inequality (4.1.3b) is equivalent to the following con-
dition.
For each connected component Γj∈Γwith Γj∩Γc6=∅there is ηj∈R
such that the following pointwise representation holds for a.a. x∈Γc
¯g(x) = PGad {−1
α(−(ν∂λ
∂n−πn)(x)−˜
F0
~el,g(¯u, ¯g)(x)
+ξ˜
F0
~ed,g(¯u, ¯g)(x) + ηjn(x))}(4.1.7)
for x∈Γc∩Γjand the zero net-mass conditions
ZΓc∩Γj
¯g·nds= 0 ∀j: Γj∩Γc6=∅
are satisfied. Here, PG:Rn→Rndenotes the Euclidean projection in Rn
onto the set G.
Proof. At first, the variational inequality (4.1.3b) is equivalent to
¯g=PGad −1
α−ν∂λ
∂n−πn−˜
F0
~el,g(¯u, ¯g) + ξ˜
F0
~ed,g(¯u, ¯g),
32
4.2. SECOND-ORDER SUFFICIENT OPTIMALITY CONDITION
where PGad :L2(Γ) →H0(Γ) is the projection with respect to the L2(Γ)-norm
on Gad. That is, ¯gsolves the minimization problem
min 1
2
g+1
α−ν∂λ
∂n−πn−˜
F0
~el,g(¯u, ¯g) + ξ˜
F0
~ed,g(¯u, ¯g)
2
L2(Γ)
subject to
ZΓc∩Γj
g·nds= 0,∀j: Γj∩Γc6=∅,
g(x)∈Gad a.e. on Γc.
Then there exist Lagrange multipliers ηjassociated to the integral constraints
in this auxiliary problem, and the variational inequality
α¯g−ν∂λ
∂n−πn−˜
F0
~el,g(¯u, ¯g) + ξ˜
F0
~ed,g(¯u, ¯g) + X
j
ηjχjn, g −¯g
≥0
(4.1.8)
holds for all g∈L2(Γ) satisfying g(x)∈Gad a.e. on Γc. By standard
arguments [102], it can be proven that this variational inequality is equivalent
to the projection representation as claimed.
Remark 4.4. The derivative of ˜
F~ei,i∈ {d, l}with respect to gis
˜
F0
~ei,g(u(x), g(x)) = −ZΓ
νg(x)∂ϕi
∂n(x)+(g(x)·n(x))(ϕi(x)·g(x)) ds.
Unfortunately, the argument −˜
F0
~el,g(¯u, ¯g) + ξ˜
F0
~ed,g(¯u, ¯g)of the projection
depends on the control itself. This term involves no smoothing operation.
Hence, we cannot conclude higher regularity of optimal controls from the
projection representation, such that ¯uhas the same regularity as ν∂λ
∂n−πn.
Moreover, the non-smooth formulation of the variational inequality is not
suitable for semi-smooth Newton methods.
To handle this difficulty, we will consider in the section after the next an
optimal control problem allowing more regular control functions. The idea
behind this is to consider a finite-dimensional control space and the once
reformulated functionals ˆ
F~ei,i∈ {d, l}, see Chapter 3.
4.2 Second-order sufficient optimality condition
The presentation of second-order sufficient optimality conditions in this sub-
section follows [21, 22, 48, 103, 109].
33
CHAPTER 4. OPTIMALITY CONDITIONS
Let (¯u, ¯p, ¯g)be a fixed admissible pair that fulfills the first-order necessary
optimality condition of Theorem 4.1 together with the adjoint state ¯
λand
the Lagrange multipliers ¯
ξ, ¯ηand let us define for simplification
j:= α¯g−ν∂¯
λ
∂n−¯πn−˜
F0
~el,g(¯u, ¯g) + ξ˜
F0
~ed,g(¯u, ¯g) + X
j
¯ηjχj!n.(4.2.1)
For ε > 0, we define by
Γε,i := {x∈Γ : |ji(x)|> ε}
the set of strongly active control constraints for ¯g.
Remark 4.5. Note that the variational inequality (4.1.8) determines the
optimal control ¯guniquely on Γε,i. We obtain
¯gi(x) = ga,i(x),if ji(x)≥ε
and
¯gi(x) = gb,i(x),if ji(x)<−ε.
Following Casas, Tröltzsch and Unger [22], based on Maurer and Zowe
[70], the linearized cone L(¯u, ¯p, ¯g)is made up of all (z, µ, h)∈U×W−1,2n/(n−1)×
L2(Γ) satisfying the following conditions (4.2.2)-(4.2.4):
−ν∆z+ (¯u·∇)z+ (z·∇)¯u+∇µ= 0 in Ω
div z= 0 in Ω
z=hon Γc,
z= 0 on Γ\Γc,
(4.2.2)
˜
F0
~ed(¯u, ¯p, ¯g)(z, µ, h)≤0,(4.2.3)
h=g−¯g, g ∈Gad.(4.2.4)
Let us denote by
Pε:L2(Γc)→L2(Γc), g 7→ χΓ\Γεg
the projection operator. That means
(Pεg)(x) = g(x)on Γ\Γε
0on Γε.
34
4.2. SECOND-ORDER SUFFICIENT OPTIMALITY CONDITION
Then, we split for all v∈L(¯v)the control function ginto g1= (g−
Pεg)and g2= (Pεg).The solutions (ui, pi), i = 1,2,of the linearized state
equations (4.2.2) are associated to gi, i = 1,2.This means
v=v1+v2= (u1, p1, g1)+(u2, p2, g2).(4.2.5)
We assume that ¯v= (¯u, ¯p, ¯g)fulfills with the Lagrange multipliers ¯
l=
(¯
λ, ¯π, ¯
ξ, ¯η)the following coercivity assumption L00(¯v, ¯
l)called (SSC):
(SSC)
There exist ε > 0and δ > 0such that
Lvv(¯v, ¯
l)[(u2, p2, g2)2]≥δkg2k2
L2(Γc)
holds for all pairs (u2, p2, g2)constructed by (4.2.5).
Theorem 4.6. Let (¯v)be an admissible non-singular point for the optimal
control problem and fulfill the first-order necessary optimality condition of
Theorem 4.11 with associated λ, ξ. Assume furthermore that the coercivity
assumption (SSC) is satisfied. Then there exist α > 0and τ > 0such that
J(v)≥J(¯v) + αkg−¯gk2
L2(Γc)
holds for all admissible pairs vwith kg−¯gkL∞(Γc)≤τ.
To prove this theorem, we establish the following two lemma.
Lemma 4.7. For all g∈Gad there holds
Z
Γc
(j) (g−¯g)dx≥εkg−¯gk2
L1(Γε).(4.2.6)
Proof. Let g∈Gad and because (¯v, ¯
λ, ¯π, ¯η, ¯
ξ)with ¯v= (¯u, ¯p, ¯g)fulfill the
first-order necessary optimality condition (4.11), we have
(ji(x)) (gi(x)−¯gi(x)) ≥0(4.2.7)
for almost all x∈Γc,i= 1, . . . , n. Integrating (4.2.7) over Γcleads with the
definition of Γε,i to
Z
Γc
(ji(x)) (gi(x)−¯gi(x))dx
≥Z
Γε,i
(ji(x)) (gi(x)−¯gi(x))dx
=Z
Γε,i
|ji(x)||gi(x)−¯gi(x)|dx
≥εZ
Γε,i
|gi(x)−¯gi(x)|dx
35
CHAPTER 4. OPTIMALITY CONDITIONS
and the sum of i= 1,2to the statement of the theorem.
Lemma 4.8. For all v=v1+v2defined as in (4.2.5) in the linearized cone
L(¯v)there holds
Lvv(¯v, ¯
l)[v1, v2]≤ckg1kL2(Γc)kg2kL2(Γc).
Proof. We get
Lvv(¯v, ¯
l)[v1, v2] = ZΩ−(u1·∇)¯
λu2)−(u2·∇)¯
λu1)dx
+ZΩ
(u1·∇)ϕlu2+ (u2·∇)ϕlu1dx
+ZΩ
(u1·∇)ϕdu2+ (u2·∇)ϕdu1dx
+ZΓ
(g1·n)(ϕl·g2)+(g2·n)(ϕl·g1)ds
+ZΓ
(g1·n)(ϕd·g2)+(g2·n)(ϕd·g1)ds+α(g1, g2)L2(Γc)
≤c(ku1kUku2kU+kg1kL2(Γc)kg2kL2(Γc)).
Because (¯u, ¯g)is non-singular and ui, i = 1,2,are the solutions of the lin-
earized Navier-Stokes equations, we obtain
Lvv(¯v, ¯
l)[v1, v2]≤ckg1kL2(Γc)kg2kL2(Γc).
Now, we are able to proof Theorem 4.6. This proof is based on [22, 109].
Proof of Theorem 4.6. We assume that (¯u, ¯p, ¯g)fulfill the assumptions of the
theorem and let (u, p, g)∈ O(¯u)×O(¯p)×O(¯g)be another admissible pair.
Then we have with v= (u, p, g)and ¯
l= (¯
λ, ¯π, ¯
ξ, ¯η)
J(¯v) = L(¯v, ¯
l)−¯
ξ(F~ed(¯v)−D0) = L(¯v, ¯
l),
J(v) = L(v, ¯
l)−¯
ξ(F~ed(v)−D0)≥ L(v, ¯
l)
due to the complementary condition. This yields
J(v)−J(¯v)≥ L(v, ¯
l)−L(¯v, ¯
l).
A Taylor-expansion of the Lagrangian Lyields the following equality:
L(v, ¯
l) = L(¯v, ¯
l) + ∂L
∂(u, p)(¯v, ¯
l)(u−¯u, p −¯p)
+∂L
∂g (¯v, ¯
l)(g−¯g) + 1
2Lvv(¯v, ¯
l)[(v−¯v)]2.
(4.2.8)
36
4.2. SECOND-ORDER SUFFICIENT OPTIMALITY CONDITION
Due to the quadratic nature of the nonlinear term, the remainder vanishes
and because the first-order necessary optimality conditions are satisfied at ¯v
with the corresponding Lagrange multipliers ¯
l, the second term of (4.2.8) is
equal to zero. For the third term we have the inequality
∂L
∂g (¯v, ¯
l)(g−¯g) = Z
Γc
(j) (g−¯g)dx≥εkg−¯gk2
L1(Γε)
see Lemma 4.7. This leads to
J(v) = J(¯v) + ∂L
∂(u, p)(¯v, ¯
l)(u−¯u, p −¯p)
+∂L
∂g (¯v, ¯
l)(g−¯g) + 1
2Lvv(¯v, ¯
λ, ¯π, ¯
ξ)[(v−¯v)]2
≥J(¯v) + εkg−¯gk2
L1(Γε)+1
2Lvv(¯v, ¯
l)[(v−¯v)]2.
Analogous to Casas, Tröltzsch and Unger [22, Section 7.2, proof of Theorem
4.2], we approximate v−¯vby vl= (ul, pl, gl)of the linearized cone L(¯v)and
the remainder term e1= (v−¯v)−vlsatisfies the estimate
ke1k ≤ ckg−¯gk2
L2(Γc).
Let us now take ul+e1instead of u−¯u, then we derive
∂2L
∂(u, p)2(¯v, ¯
l)[u−¯u, p −¯p]2=∂2L
∂(u, p)2(¯v, ¯
l)[ul, pl]2
+ 2 ∂2L
∂(u, p)2(¯v, ¯
l)[ul, pl, e1]
+∂2L
∂(u, p)2(¯v, ¯
l)[e1]2
=∂2L
∂(u, p)2(¯v, ¯
l)[ul, pl]2+e2,
with the remainder e2estimated due to the non-singularity of (¯u, ¯g)by
|e2| ≤ c(kulkU+kplkW−1,2n/(n−1) +ke1kU)ke1kU≤c(kglkL2(Γc)+ke1kU)ke1kU
and fulfilling
|e2|
kglk2
L2(Γc)→0,for kglkL2(Γc)→0.
37
CHAPTER 4. OPTIMALITY CONDITIONS
Summarized, we obtain
J(v)−J(¯v)≥1
2Lvv(¯v, ¯
l)[ul, pl, gl]2+εkg−¯gk2
L2(Γε)+e2.(4.2.9)
Now, as in (4.2.5), we split vl= (ul, pl, gl)into
ul=u1+u2, pl=p1+p2, gl=g1+g2,
where (u1, p1)and (u2, p2)are the solutions of
−ν∆u1+ (¯u·∇)u1+ (u1·∇)¯u+∇p= 0 in Ω
div u1= 0 in Ω
u1=g1on Γ.
and −ν∆u2+ (¯u·∇)u1+ (u1·∇)¯u+∇p= 0 in Ω
div u2= 0 in Ω
u2=g2on Γ.
respectively, and (u2, p2, g2)fulfill (u2, p2, g2)∈L(¯v)and (g2)i= 0 on Γε,i,
i= 1,2.Thus, (SSC) applies to Lvv(¯v, ¯
l)(u2, g2).
Considering the first term of the right-hand side of (4.2.9)
Lvv(¯v, ¯
l)[(ul, gl)]2=Lvv(¯v, ¯
l)[v1]2
+ 2Lvv(¯v, ¯
l)[v1, v2]
+Lvv(¯v, ¯
l)[v2]2,
(4.2.10)
we are able to use (SSC) and obtain
Lvv(¯v, ¯
l)[(v2)]2≥δkg2k2
L2(Γc).(4.2.11)
Lemma 4.8 leads to
Lvv(¯v, ¯
l)[v1, v2]≤ckg1kL2(Γc)kg2kL2(Γc).
We estimate Lvv(¯v, ¯
l)[v1]2analogously to Lemma 4.8 by
Lvv(¯v, ¯
l)[v1, v2]≤ckg1k2
L2(Γc).
Then, we obtain
|2Lvv(¯v, ¯
l)[v1, v2] + Lvv(¯v, ¯
l)[v2]2|
≥ −ckg1kL2(Γc)(kg1kL2(Γc)+kg2kL2(Γc))
≥ −δ
2kg2k2
L2(Γc)−ckg1k2
L2(Γc),
(4.2.12)
38
4.3. FINITE-DIMENSIONAL CONTROL SET
because of Lemma 4.8 and the Lipschitz continuity of the solution mapping
of the linearized system. Considering the relation kg2k2
L2(Γc)≥1
2kglk2
L2(Γc)−
kg1k2
L2(Γc), (4.2.10), (4.2.11) and (4.2.12), we get
Lvv(¯v, ¯
l)[(vl)]2≥δ
2kg2k2
L2(Γc)−ckg1k2
L2(Γc)≥δ
4kglk2
L2(Γc)−ckg1k2
L2(Γc)
=δ
4kg−¯gk2
L2(Γc)−ckg−¯gk2
L2(Γε).
Now, we proved
J(v)−J(¯v)≥δ
8kg−¯gk2
L2(Γc)+εkg−¯gkL1(Γε)−ckg−¯gk2
L2(Γε)+e2.
Furthermore, we obtain
J(v)−J(¯v)≥δ
8kg−¯gk2
L2(Γc)+ (ε−ckg−¯gkL∞(Γε))kg−¯gkL1(Γε)+e2.
Choosing τsufficiently small, it holds
J(v)−J(¯v)≥δ
16kg−¯gk2
L2(Γc).
4.3 Finite-dimensional control set
We will now consider a regularized version of the optimal control problem
(3.3.1)–(3.3.4). In particular, the controls will now be taken from H1/2(Γ),
which leads to higher regularity of the associated states. For H1/2(Γ)-controls
one has the following regularity result, see Theorem 2.15 or [44, Theorem
VIII.4.1].
Lemma 4.9. For every g∈H1/2(Γ) the very weak solution ubelongs to
H1(Ω), and the trace of ucoincides with the control gon the boundary Γ.
Because in this case the state uis of the space H1(Ω), we are able to use
for the boundary integrals (1.0.2) and (1.0.3) the reformulation
ˆ
F~ei(u) := −ZΩ
(ν∇u·∇ϕi+ (u·∇)yϕi)dx, i ∈ {d, l},
39
CHAPTER 4. OPTIMALITY CONDITIONS
see (3.1.2), instead of ˜
F~ei,i=d, l. Then ˆ
F~ei,i=d, l is twice continuously
Fréchet differentiable from H1(Ω) to R. And it holds
ˆ
F~ei(u) = ZΓν∂u
∂n−pnϕids=˜
F~ei(u, g), i ∈ {d, l},
for smooth states uassociated to controls g. Then, we have to redefine
G2(u, g)in (4.1.1) by
G2(u, g) = ( ˆ
F~ed(u)−D0).
.
Additionally, let us introduce a finite-dimensional control space. Let ei,
i= 1 . . . l, be linearly independent functions from H1/2(Γ) with support on
Γc. Let Qad ∈Rlbe the set of admissible coefficients
Qad := {q∈Rl:qa,i ≤qi≤qb,i}, qa, qb∈Rl.
Then, we define the set of admissible controls as
ˆ
Gad,Q := (g∈H1/2(Γ) : g=
l
X
j=1
qjej, q ∈Qad).
The idea behind this space is that in many applications only a few param-
eters can be optimized. For example, in [89] the actuation in a separation
control investigation consists of a loudspeaker, where the free optimization
parameters were frequency and amplitude. So, only two parameters appear.
Instead of this construction, we could have added a penalization term like
βkgkH1/2to the cost functional. However, this additional term is not justi-
fied physically. Moreover, the optimality system would involve a variational
inequality with a non-local differential operator on Γc.
Now, we are considering the following optimization problem, henceforth
called (Pl): Minimize
J(u, g) := −ˆ
F~el(u) + α
2kgk2
L2(Γc)
subject to the very weak form of
−ν∆u+ (u·∇)u+∇p= 0 in Ω
div u= 0 in Ω
u=gon Γc,
u= 0 on Γ\Γc,
40
4.3. FINITE-DIMENSIONAL CONTROL SET
the control constraints
g∈ˆ
Gad,Q
and the integral state constraint
ˆ
F~ed(u)≤D0.
This means
g(x) =
l
X
j=1
qjej(x), q ∈Qad
and
ˆ
F~ed(u) = −Z
Ω
(ν∇u·∇ϕi+ (u·∇)uϕi)dx≤D0.
Due to the same reasons as above, existence of solutions cannot be proven
directly. Here, we would have to work with similar modifications to (Pl)as in
Section 3.3.1 above. Rather, we would like to derive a first-order optimality
system. To this end, let us assume that (¯u, ¯g)is a non-singular and locally
optimal solution of (Pl). Moreover, let us assume that a linearized Slater
point for the state constraint exists, similarly defined as in (4.1.2)
Then one can argue as above to obtain:
Theorem 4.10. Let (¯u, ¯g)be a non-singular local optimal solution for (Pl).
Let us assume that there are ˜g∈Gad,Q ∩O(¯g)and the associated ˜u∈ O(¯u)
such that the linearized Slater condition (4.1.2) is satisfied. Then there exists
a multiplier ξ≥0, an adjoint state λ∈H1
0(Ω) ∩W2,r(Ω)n, and an adjoint
pressure π∈W1,r(Ω),r∈[2,∞), such that (λ, π)is the weak solution of
−ν∆λ+ (∇¯u)Tλ−(¯u·∇)λ+∇π= (∇¯u)T(ϕl−ξϕd)−ν∆(ϕl−ξϕd)
−(¯u·∇)(ϕl−ξϕd)in Ω
div λ= 0 in Ω
λ= 0 on Γ,
(4.3.1a)
and such that the variational inequality
α¯g−ν∂λ
∂n−πn, g −¯gL2(Γc)≥0∀g∈ˆ
Gad,Q,(4.3.1b)
and the complementarity condition
ξ(ˆ
F~ed(¯u)−D0)=0, ξ ≥0,ˆ
F~ed(¯u)≤D0(4.3.1c)
are satisfied.
41
CHAPTER 4. OPTIMALITY CONDITIONS
Here, the mapping ¯g7→ (ν∂λ
∂n−πn)is differentiable from L2(Γc)to
L∞(Γc). Thus, we are able to apply the semi smooth Newton method to
the system of Theorem 4.10. This was not possible for the system in The-
orem 4.1 due to the lack of smoothing in the argument of the projection
(4.1.7).
The variational inequality (4.3.1b) can be written equivalently as a vari-
ational inequality for the coefficients ¯qof ¯g. Let us define the mass matrix
Mand a vector Das:
M∈Rl,l, Mi,j =Z
Γc
eiejds, D ∈Rl, Di=Z
Γcν∂λ
∂n−πneids.
Then (4.3.1b) is equivalent to
(αM ¯q−D)T(q−¯q)≥0∀q∈Qad,(4.3.2)
which is the necessary and sufficient optimality condition of the quadratic
programming problem
min
q∈Qad
α
2qTMq −DTq.
Under some additional assumptions, we can simplify the system (4.3.1)
even more. Here, we will replace the functions ϕland ϕdwith differently de-
fined functions. Let us assume that there exists functions (ϕi, πi)∈H2(Ω) ×
H1(Ω),i∈ {d, l}, such that it holds
div ϕi= 0 on Ω
ϕi=eion Γw
ϕi= 0 on Γ\Γw
ν∂ϕi
∂n−πin= 0 on Γc.
(4.3.3)
The main advantage will be that that we need not calculate ϕland ϕdin the
numerical investigation. Of course, the ϕicannot be chosen as solutions of a
Stokes system, since the above conditions represent over-determined bound-
ary conditions. With this choice of auxiliary functions, all result remain
valid, since we have never used that (ϕi, πi)should be solutions of a Stokes
equation. Introducing a new adjoint state as the difference of ϕl−ξϕdand
the adjoint state given by Theorem 4.10, we obtain
Theorem 4.11. Let the assumptions of Theorem 4.10 be satisfied. Assume
there exists (ϕi, πi)∈H2(Ω)×H1(Ω),i∈ {d, l}, such that (4.3.3) is satisfied.
42
4.3. FINITE-DIMENSIONAL CONTROL SET
Then there exists a multiplier ξ≥0, an adjoint state λ∈H1(Ω) ∩W2,r(Ω),
and an adjoint pressure π∈W1,r(Ω),r∈[2,∞), such that (λ, π)is the weak
solution of
−ν∆λ+ (∇¯u)Tλ−(¯u·∇)λ+∇π= 0 in Ω
div λ= 0 in Ω
λ=~el−ξ~edon Γw
λ= 0 on Γ\Γw,
(4.3.4a)
and such that the variational inequality
α¯g−ν∂λ
∂n−πn+ (¯u·n)λ, g −¯gL2(Γc)≥0∀g∈ˆ
Gad,Q,(4.3.4b)
and the complementarity condition
ξ(ˆ
F~ed(¯u)−D0)=0, ξ ≥0,ˆ
F~ed(¯u)≤D0(4.3.4c)
are satisfied.
The additional term (¯u·n)λin (4.3.4b) appears because of λ∈H1(Ω) ∩
W2,r(Ω) instead of H1
0(Ω) ∩W2,r(Ω),r∈[2,∞), see e.g. Section 5.2.1 for a
formal Lagrange approach.
Please note that this system does not involve the functions ϕiin the
adjoint equation and in the variational inequality, which was the case for the
optimality systems obtained in Theorems 4.1, 4.2, and 4.10. This makes this
system favorable for computations, and it is used for the solution algorithm
that we employed in our numerical experiments.
The optimality system (4.3.4a)-(4.3.4c) of the previous theorem can also
be obtained formally with the help of the Lagrangian. Let us define the La-
grange functional Las the sum of the cost functional, Navier-Stokes equation
tested with (λ, π), and state constraint
L(u, p, g, λ, π, ξ) = −Z
Γwν∂u
∂nw−pnw·~elds+α
2Z
Γc
|g|2ds
−Z
Ω
λ(−ν∆u+ (u·∇)u+∇p)dx
+Z
Γc
(g−u)λ2ds+Z
Γc
uλ3ds
+ξ
−Z
Γwν∂u
∂n−pn·~edds−D0
.
43
CHAPTER 4. OPTIMALITY CONDITIONS
Then, (4.3.4a) is equivalent to
Lu(¯u, ¯p, ¯g, λ, π, ξ)v= 0,∀v∈H1(Ω),
Lp(¯u, ¯p, ¯g, λ, π, ξ)q= 0,∀q∈L2(Ω)
and (4.3.4b) can be obtained formally by
Lg(¯u, ¯p, ¯g, λ, π, ξ)(g−¯g)≥0,∀g∈ˆ
Gad,Q.
In order to obtain (4.3.4b) in a rigorous way, we had to use regular controls
and to suppose the existence of ϕisatisfying (4.3.3).
We have chosen
ˆ
Gad,Q := (g∈H1/2(Γ) : g=
l
X
j=1
qjej, q ∈Qad)
with ei∈H1/2(Γ). Let us consider for simplicity in the following chapters
Gad,Q := (g∈H3/2(Γ) : g=
l
X
j=1
qjej, q ∈Qad)(4.3.5)
with ei∈H3/2(Γ) and qa, qb∈Rl.
Based on this set of admissible controls Gad,Q, we get u∈H2(Ω), see
Theorem 2.17.
Depending on the respective situation in the following, we will decide
to apply one of the two equivalent notations (4.3.2) or (4.3.1b) with Gad,Q
instead of ˆ
Gad,Q.
Now, we are also able to use F~ei(u)in (4.3.4c) and the cost functional
instead of ˆ
F~ei(u),i∈ {d, l}.Additionally, we have to redefine G2(u, g)in
(4.1.1) by
G2(u, g) = (F~ed(u)−D0).
.
Furthermore, we have to redefine the definition of non-singularity and
Theorem 2.26 with (¯u, ¯g)∈H2(Ω)×H3/2(Γc)or (¯u, ¯q)∈H2(Ω)×Rl, respec-
tively.
4.4 Second-order sufficient optimality conditions
for the finite-dimensional case
Almost analogously to the infinite-dimensional case, we are able to consider
the second-order sufficient optimality conditions with the control space Qad.
We refer also to [31] for this case.
44
4.4. SECOND-ORDER SUFFICIENT OPTIMALITY CONDITIONS
FOR THE FINITE-DIMENSIONAL CASE
Let us therefore define d:= (αM ¯q−D)T, similar to (4.3.2), with the mass
matrix Mand a vector D
M∈Rl,l, Mi,j =Z
Γc
eiejds, D ∈Rl, Di=Z
Γcν∂λ
∂n−πn+ (¯u·n)λeids.
Then (4.3.4b) is equivalent to
(αM ¯q−D)T(q−¯q)≥0∀q∈Qad,(4.4.1)
which is the necessary and sufficient optimality condition of the quadratic
programming problem
min
q∈Qad
α
2qTMq −DTq.
Additionally, we introduce A+:= {i:di>0},A−:= {i:di<0}and
A:= A+∪A−and the critical cone associated with ¯q
C¯q:= {h∈Rl:hi= 0 ∀i∈ A satisfying (4.4.2) −(4.4.4)}
−ν∆z+ (¯u·∇)z+ (z·∇)¯u+∇µ= 0 in Ω
div z= 0 in Ω
z(x) =
l
X
j=1
ej(x)hjon Γc,
z= 0 on Γ\Γc,
(4.4.2)
(F0
~ed(¯u)(z) = 0,if F~ed(¯u) = D0and ξ > 0
F0
~ed(¯u)(z)≤0,if F~ed(¯u) = D0and ξ= 0,,(4.4.3)
hi=(≥0if ¯qi=qa,i
≤0if ¯qi=qb,i
.(4.4.4)
Then, we define the coercivity assumption (SSC0):
(SSC0)
The inequality
hTLqq(¯q, ¯
ξ)h > 0
holds for all h∈C¯q\{0}.
and derive the following theorem
Theorem 4.12. Let (¯u, ¯g)be an admissible non-singular point for the opti-
mal control problem and fulfill the first-order necessary optimality condition
of Theorem 4.11 including the variational inequality (4.3.2) with associated
45
CHAPTER 4. OPTIMALITY CONDITIONS
λ, ξ. Assume furthermore that the coercivity assumption (SSC’) is satisfied.
Then there exist δ > 0and τ > 0such that
J(v)≥J(¯v) + δ|q−¯q|2
holds for all admissible pairs (u, g)with |q−¯q| ≤ τ.
46
Chapter 5
Numerical investigations
In this chapter, we want to provide numerical results for the problem under
consideration. We handle the optimal control problem above numerically by
direct solution of the optimality system that follows from Theorem 4.11 and
is stated below. We follow a method suggested by Neitzel et al. [73].
Afterwards, we consider an SQP-method, whose convergence is proved in
Chapter 6.
5.1 One-shot approach
Here, we consider a slightly different problem. We have an inflow g∞acting
as an inhomogeneous Dirichlet boundary condition on the inflow boundary
Γin. The control boundary Γcwas modelled by a nonhomogeneous Dirichlet
condition, where the limited suction and/or blowing occurs on small slot on
the flap. A no-slip boundary condition, i.e. homogeneous Dirichlet condition,
was used for the remaining airfoil Γwand the wall boundaries Γwall. At the
outflow Γout, we applied a so called ’do-nothing’-condition:
ν∂u
∂n−pn= 0.
For more details of the configuration see the technical report [19].
These do-nothing conditions have similar properties as Neumann bound-
ary conditions for scalar elliptic equations. Let us briefly comment on avail-
able results for Navier-Stokes equations with mixed boundary conditions of
Dirichlet and Neumann type. For the Navier-Stokes system with mixed
boundary conditions, existence and uniqueness of solutions u∈H1(Ω) for
Dirichlet data g∈H1/2(Γc)were proven in [72], Theorem 5.2, and [71]. We
can show, similar to Theorem 2.17, that we obatin a solution u∈H2(Ω) for
47
CHAPTER 5. NUMERICAL INVESTIGATIONS
g∈H3/2(Γc). Due to the ’do-nothing’ boundary condition on the outflow
boundary, the pressure p∈L2(Ω) is unique. To the best of our knowledge,
similar results for low-regularity Dirichlet data in L2(Γc)are still missing. If
one carries out the formal procedure as described at the end of the finite-
dimensional part or see (5.2.5), one finds that the adjoint equation for the
problem with ’do-nothing’ outflow condition is given by
−ν∆λ+ (∇¯u)Tλ−(¯u·∇)λ+∇π= 0 in Ω
div u= 0 in Ω
λ=~el−ξ~edon Γw
λ= 0 on Γ\(Γw∪Γout)
ν∂λ
∂n−πn+ (¯u·n)λ= 0 on Γout.
This adjoint system is analogously to the one of Theorem 4.11.
As already mentioned above, we solved the optimality system given anal-
ogously to Theorem 4.11. Due to the presence of the ’do-nothing’ boundary
condition, we can drop the constraint RΓjg·nds= 0, which was incorporated
to guarantee existence of divergence-free solutions. With this simplification,
the variational inequality and the complementarity condition in the optimal-
ity system given by Theorem 4.11 are equivalent to
¯g=PG1
α(ν∂λ
∂n −πn + (¯u·n)λ)a.e. in Γc
and
ξ= max (0, ξ +F~ed(¯u)−D0).
This enables us to eliminate the control variable by means of the projection.
Then we want to solve the following system consisting of the state equa-
tion with control eliminated by the projection formula, the equation for the
48
5.1. ONE-SHOT APPROACH
Lagrange multiplier ξand the associated adjoint equation, see also [19]:
−ν∆u+ (u·∇)u+∇p= 0 in Ω
div u= 0 in Ω
u=PG{1
α(ν∂λ
∂n −πn + (u·n)λ)}on Γc
u= 0 on Γwall ∪Γw\Γc
u=g∞on Γin
ν∂u
∂n −pn = 0 on Γout
−ν∆λ+ (∇u)Tλ−(u·∇)λ+∇π= 0 in Ω
div λ= 0 in Ω
λ=~el−ξ~edon Γw
λ= 0 on Γin ∪Γwall
ν∂λ
∂n −πn + (u·n)λ= 0 on Γout
and
ξ= max (0, ξ +F~ed(u)−D0).(5.1.1)
5.1.1 Numerical results
The computational domain, depicted in Figure 5.1, is a 2D generic high-lift
configuration consists of a NACA4412 main airfoil at 8◦angle of attack and
a NACA4415 flap with a deflection angle of 37◦. The Reynolds number was
given as Re = 85 based on the chord length Lref = 1.275 and the free stream
velocity g∞= 1.
We used the commercial FEM-solver COMSOL Multiphysics with a build-
in damped Newton method for the nonlinear system. The partial differential
equations were discretized using Taylor-Hood finite elements, i.e. piecewise
quadratic polynomials for the velocity and piecewise linear polynomials for
the pressure.
The equation for ξis in this form not solvable, because on the one hand
we only can define variables on the whole area Ωand on the other hand
COMSOL does not allow variables in R.Therefore, we choose the following
algorithm. In the first step, we solve the optimal control problem without the
state constraint. If for the computed solution the state constraint is satisfied
then this solution solves also the state-constrained problem. Otherwise, we
49
CHAPTER 5. NUMERICAL INVESTIGATIONS
Figure 5.1: The domain
have to consider the state constraint
Z
Γw
ν∂u
∂n−pn·~edds=D0.
Integral terms in the system and COMSOL Multiphysics don’t fit together.
To handle this problem, we consider an augmented Lagrangian method for
this problem, see [13, 46, 78]. This algorithm is related to the Penalty prob-
lem. The difference is that we here reduce the risk of ill-conditioned sub-
problems, because now we introduce Lagrange multiplier estimates at each
step to the cost functional. The penalty term does not guarantee
F~ed(u, p) = D0.(5.1.2)
It only leads to
(F~ed(¯u, ¯p)−D0) = −1
cξ, (5.1.3)
see [78], (17.45). So, we see (5.1.2) is theoretical fulfilled for c→ ∞, but we
probably get ill conditioned and numerical problems for big values of c.
Let us ignore for a while the Navier-Stokes equations and just consider
the cost functional, which we want to minimize, subject to the integral state
50
5.1. ONE-SHOT APPROACH
constraint. Then the problem is
(P1) min
g∈Gad,Q
J(u, p, g)subject to F~ed(u, p) = D0
with the optimality system
∇J(¯u, ¯p) + ξF0
~ed(¯u, ¯p)=0.(5.1.4)
The associated Penalty problem reads as
(P2) min
g∈Gad,q
J(u, p, g) + c(F~ed(u, p)−D0)2
with the optimality system
∇J(¯u, ¯p)+2cF0
~ed(¯u, ¯p)(F~ed(¯u, ¯p)−D0)=0 (5.1.5)
The augmented Lagrangian function LA(u, g, ξ)avoids the problem that
we need very big penalty parameter cby an estimation for the Lagrange-
multiplier ξfor the integral state constraint. The augment Lagrangian
LA(u, g, ξ)consists of the original cost functional J(u, g), the penalty term
and the term involving the multiplier ξ:
LA(u, g, ξ) := J(u, g)−ξ(F~ed(u, p)−D0) + c(F~ed(u, p)−D0)2.
Considering LA(u, g, ξ)as the new cost functional, we obtain the optimality
system
∇(u,p)J(u, g)−F0
~ed(¯u, ¯p)(ξ−2c(F~ed(¯u, ¯p)−D0)) = 0.(5.1.6)
Now, we fix the penalty parameter cand the Lagrangian multiplier ξin
each step by cand ξk, respectively. The optimality systems (5.1.4), (5.1.5)
and (5.1.6) leads to
ξ≈ξk−2c(F~ed(¯u, ¯p)−D0)(5.1.7)
which is equivalent to
(F~ed(¯u, ¯p)−D0)≈1
2c(ξk−ξ)(5.1.8)
and we see that if ξkis close to the original Lagrange multiplier ξ, (5.1.8) is
closer to F~ed(¯u, ¯p) = D0than (5.1.3). From (5.1.7), we get also the prescript
to calculate the new Lagrange multiplier in the k+ 1th step by
ξk+1 =ξk−2c(F~ed(¯u, ¯p)−D0).
51
CHAPTER 5. NUMERICAL INVESTIGATIONS
Arada, Raymond and Tröltzsch proved in [13] the convergence of the aug-
mented Lagrangian method for a class of problems.
If the integral state constraint is fulfilled, then we are able to reduce the
associated Lagrange-multiplier ξ. Or otherwise, we have to increase the mul-
tiplier.
For the uncontrolled problem, we obtained a lift of CA=FA
0.5g2
∞Lref = 1.562
and a drag of CD=DA
0.5g2
∞Lref = 0.817, where FAis the resulting lift force,
DAthe drag force and Lref = 1.275 is the reference length of the wing, see
Figure 5.2 for a streamline plot of the velocity field, and Figure 5.4 for a plot
of the absolute values of the uncontrolled velocity field and Figure 5.5 for the
pressure field.
Figure 5.2: Uncontrolled case: velocity field.
52
5.1. ONE-SHOT APPROACH
Figure 5.3: Uncontrolled case: velocity field with zoom on the wing.
Figure 5.4: Uncontrolled case: absolute value of velocity field (left).
53
CHAPTER 5. NUMERICAL INVESTIGATIONS
Figure 5.5: Uncontrolled case: absolute value of pressure field (right).
Now let us report about the result of the optimization. Here, we choose
the control cost parameter α= 0.1and the control constraints as box con-
straints G= [−1,+1].
At first, we compute the solution for the case without any drag constraint.
The optimal control is given by the maximal possible suction, which is natural
from a physical point of view. The obtained optimized lift is CA= 1.5823 and
the drag is CD= 0.8340, which is a lift gain of 1.3%. The controlled velocity
field can be seen in Figure 5.6. The adjoint velocity field and pressure are
plotted in Figures 5.8 and 5.9.
54
5.1. ONE-SHOT APPROACH
Figure 5.6: Controlled case: velocity field.
Figure 5.7: Controlled case: velocity field with zoom on the wing.
55
CHAPTER 5. NUMERICAL INVESTIGATIONS
Figure 5.8: Controlled case: adjoint velocity field
Figure 5.9: Controlled case: absolute value of adjoint velocity field (left).
56
5.2. SQP-METHOD
Figure 5.10: Controlled case: image of the slit.
The box constraint G= [−0.5,+0.5] leads to a lift coefficient of CA=
1.572 and a drag coefficient of CD= 0.825, which amounts to a lift gain of
about 0.65%.
In the next step, we choose D0= 0.5240 as upper bound for the drag.
This equates to an constraint coefficient of CD0= 0.8220. Hence, we expect
that this constraint will be active at the solution. In fact, for the com-
puted solution we obtain CD= 0.8215. Moreover, due to this restriction the
computed lift is CA= 1.571, which is less than for the case without state
constraints, but which is still better than in the uncontrolled situation. An
upper bound for the drag of D0= 0.5230,CD0= 0.820 leads to Cd= 0.820.
and Ca= 1.570.
An upper bound for the drag of D0= 0.5190,CD0= 0.814 leads to
Cd= 0.814. and Ca= 1.559.
5.2 SQP-method
In addition to the last section, where we solved the optimality system at
once, we will now solve the problem by an iterative method. A widely-used
method is the SQP-method.
57
CHAPTER 5. NUMERICAL INVESTIGATIONS
5.2.1 The SQP-method for our problem
In this case, we do not solve the original problem. The idea is to solve a
slightly different problem PI, where the integral state constraint is added as
penalty term in the cost functional. We do this, because we want to avoid
the integral equation for the Lagrange multiplier ξ, see (5.1.1).
Analogously to the one-shot approach, we first solve the problem without
a state constraint, this means in this case without the penalty term, and if
the state constraint is satisfied for the computed solution then this solution
also solves the state-constraint problem. Otherwise, it is sufficient to consider
the state constraint as an equality
ZΓwν∂u
∂n−pn·~edds=D0.
So, we can avoid the penalty term c R
Γwν∂u
∂n−pn·~edds−D0!2
+
, where
(s)+:= max{0, s},s∈R, denotes the positive part. The problem is that
this term is not twice continuously Fréchet-differentiable.
In this case, we introduce the penalty term to the cost functional and
regard the following problem: Minimize
JI(u, g) := Z
Γwν∂u
∂n−pn·~elds
+c
Z
Γwν∂u
∂n−pn·~edds−D0
2
+α
2kgk2
L2(Γc)
(5.2.1)
subject to
−ν∆u+ (u·∇)u+∇p= 0 in Ω
div u= 0 in Ω
u=gon Γc
u= 0 on Γwall ∪Γw\Γc
u=g∞on Γin
ν∂u
∂n−pn= 0 on Γout
(5.2.2)
the control constraints
g∈Gad,Q (5.2.3)
58
5.2. SQP-METHOD
where c∈Ris the penalty constant and g∞∈Rthe inflow. The solution
uis defined in a (very) weak sense. We have shown in Section 2, Theorem
2.17, that we get for every g∈H3/2(Γc)a solution u∈H2(Ω) for the state
equation (5.2.2). The optimality system could be derived very similar to the
problem with the integral state constraint. The Lagrangian looks as follows:
L(u, p, g, λ) = Z
Γwν∂u
∂n −pn·~elds+α
2Z
Γc
|g|2ds
+c
Z
Γwν∂u
∂n −pn·~edds−D0
2
+Z
Ω
(−ν∆u+ (u·∇)u)λ+∇pλdx−Z
Ω
πdiv udx+Z
Γc
(u−g)λ2ds
+Z
Γwall
uλ3ds+Z
Γin
(u−g∞)λ4ds+Z
Γout ν∂u
∂n −pnλ5ds.
with
Z
Ω
(u·∇)u·λdx=−Z
Ω
(u·∇)λ·u+u·λdiv udx+ZΓ
u·nu·λdγ.
The necessary condition ∂L
∂(u,p)(¯u, ¯p)(u, p)=0leads to
∂L
∂(u, p)(¯u, ¯p)(u, p) = Z
Γwν∂u
∂n−pn·~elds
+ 2c
Z
Γwν∂¯u
∂n−¯pn·~edds−D0
| {z }
:=K(¯u,¯p))
Z
Γwν∂u
∂n−pn·~edds
+Z
Ω
−ν∆uλ −(¯u·∇)λ·u−(u·∇)λ·¯u−¯u·λdiv u+∇pλdx
(5.2.4)
−Z
Ω
πdiv udx+Z
Γc
uλ2ds+Z
Γwall
uλ3ds+Z
Γ∞
uλ4ds
+Z
Γout ν∂u
∂n −pnλ5ds+Z
Γc∪Γout
¯u·nu·λ+u·n¯u·λds= 0.
59
CHAPTER 5. NUMERICAL INVESTIGATIONS
This is with partial integration and
K(¯u, ¯p)=2c
Z
Γwν∂¯u
∂n−¯pn·~edds−D0
equivalent to
Z
Γwν∂u
∂n−pn·~elds+K(¯u, ¯p)Z
Γwν∂u
∂n−pn·~edds
−Z
Ω
(ν∇u·∇λ−(¯u·∇)λ·u−(u·∇)λ·¯u
+ ¯u·λdiv u−pdiv λ+u·∇π)dx
+Z
Γc−ν∂u
∂nλ+pnλ−πnu+unλ2ds
+Z
Γwall −ν∂u
∂nλ+pnλ−πnu+unλ3ds
+Z
Γin −ν∂u
∂nλ+pnλ−πnu+unλ4ds
+Z
Γout −ν∂u
∂n(λ−λ5) + pn(λ−λ5) + πnuds
+Z
Γc∪Γout
¯u·nu·λ+u·n¯u·λds= 0.
The equality
−Z
Ω
((u·∇)λ·¯u+ ¯u·λdiv u)dx+Z
Γc∪Γout
(u·n)(¯u·λ)ds
=Z
Ω
(∇¯u)Tλ·udx
60
5.2. SQP-METHOD
and another partial integration lead to
Z
Γwν∂u
∂n−pn·~elds+K(¯u, ¯p)Z
Γwν∂u
∂n−pn·~edds
−Z
Ω−ν∆λ+ (¯u·∇)λ+ (∇¯u)Tλ+∇πudx−Z
Ω
pdiv λdx
−Z
Γcν∂λ
∂nu−ν∂u
∂nλ+pnλ−πnu+unλ2+ ¯u·nu·λds
−Z
Γwall ν∂λ
∂nu−ν∂u
∂nλ+pnλ−πnu+unλ3ds
−Z
Γin ν∂λ
∂nu−ν∂u
∂nλ+pnλ−πnu+unλ4ds
−Z
Γout ν∂λ
∂nu−ν∂u
∂n(λ−λ5) + pn(λ−λ5) + πnu+ ¯u·nu·λds
= 0.
Taking (u, p)∈H2
0(Ω) ×H1(Ω), we obtain
Z
Ω−ν∆λ+ (¯u·∇)λ+ (∇¯u)Tλ+∇πudx−Z
Ω
pdiv λdx
+Z
Γw
(~el+K(¯u, ¯p)~ed−λ)(ν∂u
∂n−pn)ds+Z
Γwall
λ(ν∂u
∂n−pn)ds
+Z
Γin
λ(ν∂u
∂n−pn)ds+Z
Γout
(λ−λ5)(ν∂u
∂n−pn)ds= 0
so that λand πare the weak solutions of
−ν∆λ+ (∇¯u)Tλ−(¯u·∇)λ+∇π= 0 in Ω
div λ= 0 in Ω
λ=~el+K(¯u, ¯p)~edon Γw
λ= 0 on Γin ∪Γwall
λ=λ5on Γout.
(5.2.5)
61
CHAPTER 5. NUMERICAL INVESTIGATIONS
Taking (u, p)∈H2(Ω) ×H1(Ω), we obtain
Z
Γwν∂λ
∂nu−πnu+unλ2+ ¯u·nu·λds= 0,
Z
Γwall ν∂λ
∂nu−πnu+unλ3ds= 0,
Z
Γin ν∂λ
∂nu−πnu+unλ4ds= 0,
Z
Γout ν∂λ
∂nu−πnu+ ¯u·nu·λds= 0,
which implies
λ2n=−(ν∂λ
∂n−πn+ (¯u·n)λ)on Γw,(5.2.6)
λ3n=−(ν∂λ
∂n−πn)on Γwall,
λ4n=−(ν∂λ
∂n−πn)on Γin,
0 = −(ν∂λ
∂n−πn+ (¯u·n)λ)on Γout.
Now, we are able to substitute the last equation of (5.2.5) with ν∂λ
∂n+πn+
ν∂λ
∂n+πn+ ¯u·nu·λ= 0 on Γout.
The other necessary condition ∂L
∂g (¯g)(g−¯g)≥0for all g∈Gad,Q leads to
∂L
∂g (¯g)(g−¯g) = Z
Γc
α¯g(g−¯g)−(g−¯g)λ2ds≥0∀g∈Gad,Q
which with (5.2.6) is equivalent to
Z
Γcα¯g−(ν∂λ
∂n−πn+ (¯u·n)λ)(g−¯g)ds≥0∀g∈Gad,Q.(5.2.7)
Altogether, we derive the adjoint system
−ν∆λ+ (∇¯u)Tλ−(¯u·∇)λ+∇π= 0 in Ω
div λ= 0 in Ω
λ=~el+K(¯u, ¯p)~edon Γw
λ= 0 on Γin ∪Γwall
ν∂λ
∂n+πn+ ¯u·nu·λ= 0 on Γout
(5.2.8)
62
5.2. SQP-METHOD
with K(¯u, ¯p) = 2c R
Γwν∂¯u
∂n−¯pn·~edds−D0!and the variational inequal-
ity
Z
Γcα¯g−(ν∂λ
∂n−πn+ (¯u·n)λ)(g−¯g)ds≥0∀g∈Gad,Q.(5.2.9)
Let (¯u, ¯g, ¯
λ)satisfy the optimality system, consisting of the Navier-Stokes
equations (5.2.2), the control constraint (5.2.3), the adjoint system (5.2.8)
and the variational inequality (5.2.9).
First, we consider the problem without any restrictions to the control
function gthat means Gad,Q ={v∈H3/2(Γc) : v(x) = Pl
i=1 ei(x)qi, q ∈Rl}
or Qad =Rland instead of the variational inequality (5.2.9), we have the
equation α¯g−(ν∂λ
∂n+πn+ (¯u·n)λ) = 0. So, the optimality systems reads as
(5.2.2),(5.2.3),(5.2.8) and
α¯g−ν∂λ
∂n−πn+ (¯u·n)λ= 0 ∀g∈Gad,Q.
This nonlinear system for (u, g, λ)can be solved with the Newton-method,
see for instance [32]. Let us briefly sketch this method: consider the opti-
mization problem
min f(u)u∈Rn
with f∈C2(Rn).
Then, we want to solve the optimality system
f0(¯u) = 0.
Solving this system with the Newton-method means that we get un+1 as the
solution of
f0(un) + f00(un)(u−un)=0
When we transfer this idea to our problem, we obtain the following Newton-
method.
Algorithm 5.1 (NM).
1. Choose an initial value z0= (u0, g0, λ0)and set k= 0.
63
CHAPTER 5. NUMERICAL INVESTIGATIONS
2. Determine z= (u, g, λ)by the Navier-Stokes system linearized at zk
−ν∆u+ (u·∇)uk+ (uk·∇)u+∇p=−(uk·∇)ukin Ω
div u= 0 in Ω
u=gon Γc
u= 0 on Γwall ∪Γw\Γc
u=g∞on Γin
ν∂u
∂n−pn= 0 on Γout,
(5.2.10a)
−ν∆λ+ (∇¯u)Tλ−(¯u·∇)λ+∇π=−(∇(¯u−uk))Tλk
+ ((¯u−uk)·∇)λkin Ω
div λ= 0 in Ω
λ=~el+K(¯u, ¯p)~edon Γw
λ= 0 on Γin ∪Γwall
ν∂λ
∂n −πn + (¯u·n)λ= 0 on Γout.
(5.2.10b)
and
α¯g−ν∂λ
∂n−πn+ (¯u·n)λ= 0 (5.2.10c)
3. Set k=k+ 1 and zk=z. Goto 2.
The equations (5.2.10) are equivalent to solving the quadratic problem
min Jk(u, g) := ∇J(uk, gk)(u−uk, g −gk)
+1
2Lzz(uk, gk, λk)[(u−uk, g −gk)]2(5.2.11)
subject to the linearized Navier-Stokes equation
−ν∆u+ (u·∇)uk+ (uk·∇)u+∇p=−(uk·∇)ukin Ω
div u= 0 in Ω
u=gon Γc
u= 0 on Γwall ∪Γw\Γc
u=g∞on Γin
ν∂u
∂n−pn= 0 on Γout.
(5.2.12)
So, we can also solve this system instead of (5.2.10). In contrast to the
system above, the system (5.2.10) cannot be transformed to the case with
64
5.2. SQP-METHOD
g∈Gad,Q.
Instead, we have to solve in the n-th step (QPn):
min Jk(u, g) :=∇J(uk, gk)(u−uk, g −gk)
+1
2Lzz(uk, gk, λk)[(u−uk, g −gk)]2
−ν∆u+ (u·∇)uk−(uk·∇)u+∇p=−(uk·∇)ukin Ω
div u= 0 in Ω
u=gon Γc
u= 0 on Γwall ∪Γw\Γc
u=g∞on Γin
ν∂u
∂n−pn= 0 on Γout
g∈Gad,Q.
The cost functional, solved in the k-th step, Jkdiffers only by the term
1
2Lzz(uk, gk, λk)[(u−uk, g −gk)]2=(((u−uk)·∇)(u−uk), λk)L2(Ω)
+ 2c(ν∂(u−uk)
∂n−(p−pk)n,~el)L2(Ω)
from the original cost functional J. Our idea is to solve the linear subprob-
lems (QPn)with the gradient-projection method.
Remark 5.2. We know that it is not the best method; we solve the whole
problem with the (fast) SQP-method and the subproblems with the (slow)
gradient-projection method. A more appropriate choice would be for instance
the active-set method. In this case, we would have to handle the integral term
K(u, p)=2cRΓwν∂u
∂n−pn·~edds−D0. But, we decided to solve our
problems numerically with COMSOL Multiphysics. As mentioned before, this
does not fit together. Therefore, we consider the gradient-projection method
to calculate K(¯u, ¯p)after solving the state equation.
5.2.2 Gradient-projection method
Let us briefly formulate the principle of this method for an optimization
problem in a Hilbert space U,
min
u∈Uad
f(u),
where Uad ⊂Uis a non-empty, bounded, convex and closed set and f:U→
Ris a Gâteaux-differentiable functional.
65
CHAPTER 5. NUMERICAL INVESTIGATIONS
The iteration steps u1, ..., unare finished so that unis the current solution.
Then the algorithm look as follows:
S1 (Direction search) Choose the anti-gradient as descent-direction
vn:= −f0(un).
S2 (Stepsize) Choose a step size sn, so that the equation
f(P[ua,ub]{un+snvn}) = min
s>0f(P[ua,ub]{un+svn})
is fulfilled, which guarantees the admissibility of the solution.
S3 Set un+1 =un+snvn,n=n+ 1 and goto S1.
For the optimization subproblems (QPn)the algorithm reads as follows
S1 Calculate (un, pn)as the solution of (5.2.10a).
S2 Calculate the adjoint (λn, πn)from (5.2.10b).
S3 The updated descent direction is
vn:= αgn−(ν∂λ
∂n+πn+ (¯u·n)λ)
S4 Calculate the step size snfrom
min
s>0f(P[ga(x),gb(x)]{gn+svn}).
S5 The updated control gn+1 is
gn+1 := P[ga(x),gb(x)]{gn+snvn}.
set n:=n+1 and goto S1.
5.2.3 Example
Let us now consider the same example as for the one-shot approach with
ga(x)≡ −0.5,gb(x)≡0.5and a Reynolds number Re =1
νLref = 85.
With a drag constraint of D0= 0.5240 which is equal to a coefficient
CD0= 0.822 and a penalty parameter of c= 50, we obtain CA= 1.571 and
CD= 0.8215 . The results are presented in the Figures 5.11 and 5.12.
66
5.2. SQP-METHOD
An upper constraint of D0= 0.5190,CD0= 0.814 for the drag coefficient
leads to CD= 0.814 and CA= 1.559.
We recognize that we get the same results for the SQP-method than for
the one-shot approach.
Figure 5.11: Streamlines for the optimal velocity field with an integral state
constraint D0= 0.5240.
67
CHAPTER 5. NUMERICAL INVESTIGATIONS
Figure 5.12: Streamlines for the optimal velocity field with an integral state
constraint D0= 0.5240 an a zoom on the slit.
68
Chapter 6
Convergence of the SQP-method
In this chapter, we want to prove the convergence of the SQP-method, men-
tioned in Section 5.2. Our approach is mainly based on the theses of A.
Unger [104] and D. Wachsmuth [109]. Additionally, we refer to [84].
We want to prove the convergence of the SQP-method for the following
problem with penalty term: Minimize
JI(u, q) := Z
Γwν∂u
∂n−pn·~elds+α
2kgk2
L2(Γc)
+c
Z
Γwν∂u
∂n−pn·~edds−D0
2
subject to the very weak form of the nonhomogenous Navier-Stokes equations
−ν∆u+ (u·∇)u+∇p= 0 in Ω
div u= 0 in Ω
u=
l
X
i=1
eiqion Γc
u= 0 on Γ\Γc
and the control constraints
q∈Qad,
where c∈Ris the penalty constant.
The associated optimality system consists analogously to (5.2.2), (5.2.8)
69
CHAPTER 6. CONVERGENCE OF THE SQP-METHOD
and (5.2.9) of the state equations
−ν∆u+ (u·∇)u+∇p= 0 in Ω
div u= 0 in Ω
u=g=
l
X
i=1
eiqion Γc
u= 0 on Γ\Γc
the adjoint system
−ν∆λ+ (∇¯u)Tλ−(¯u·∇)λ+∇π= 0 in Ω
div λ= 0 in Ω
λ=~el+K(¯u, ¯p)~edon Γw
λ= 0 on Γ\Γw
and the variational inequality
Z
Γcα¯g−(ν∂λ
∂n−πn+ (¯u·n)λ)(g−¯g)ds≥0∀g∈Gad,Q
or equivalent
(αM ¯q−D)T(q−¯q)ds≥0∀q∈Qad
with
M∈Rl,l, Mi,j =Z
Γc
eiejds, D ∈Rl, Di=Z
Γcν∂λ
∂n−πn+ (¯u·n)λeids.
Consider with d:= (αM ¯q−D)T, see (4.4.1), the sets A+:= {i:di>0},
A−:= {i:di<0},A:= A+∪A−and the critical cone associated with ¯q
C¯q:= {h∈Rl:hi= 0 ∀i∈ A satisfying (6.0.2)}
hi=(≥0if ¯qi=qa,i
≤0if ¯qi=qb,i
.(6.0.2)
Then, we define the coercivity assumption (SSC00):
(SSC00)
The inequality
hTLqq(¯q)h > 0
holds for all h∈C¯q\{0}.
70
6.1. GENERALIZED EQUATIONS
There are several problems in handling the optimality system, especially
the nonlinearity of the state equation and the control. The Newton-method
is very popular to solve nonlinear systems of the form
0 = f(x),
because of the fast local convergence. Due to the variational inequality,
we are not able to use it in the classical form. A loophole are generalized
equations. Therefore, we have first to introduce generalized equations.
6.1 Generalized equations
Let the normal cone given by
NC(u) := {z∈Rn:zT(v−u)≤0∀v∈C},
and ˜
Gc:H3/2(Γc)→H2(Ω), g 7→ u,
Gw:H3/2(Γw)→H2(Ω), g 7→ u
and
S:L2(Ω) →H2(Ω), f 7→ u
denote the control-to-solution operators of the Stokes equations, see Corollary
2.12. Let us furthermore define the operators
Hc:Rl→H3/2(Γc), q 7→
l
X
i=1
ei(x)qi
and
Gc:= ˜
Gc◦Hc:Rl→H2(Ω), q 7→ u.
Let us recall
K(¯u, ¯p)=2c
Z
Γwν∂¯u
∂n−¯pn·~edds−D0
and we denote the derivative by
˜
K(u, p) := ∂
∂(u, p)K(¯u, ¯p)(u, p) = 2c
Z
Γwν∂u
∂n−pn·~edds
.
71
CHAPTER 6. CONVERGENCE OF THE SQP-METHOD
Then, we reformulate the optimality system (6.0.1)-(6.0.1) at ¯z= (¯u, ¯p, ¯q, ¯
λ, ¯π)
to −¯u+S(−(¯u·∇)¯u) + Gc(¯q)=0
−¯
λ+S(−(∇¯u)Tλ+ (¯u·∇)λ) + Gw(~el+K(¯u, ¯p)~ed)=0
−(αM ¯q−D)T∈ NQad (¯q)
and the equivalent formulation
0∈ {−¯u+S(−(¯u·∇)¯u) + Gc(¯q)}+{0}
0∈ {−¯
λ+S(−(∇¯u)Tλ+ (¯u·∇)λ) + Gw(~el+K(¯u, ¯p)~ed)}+{0}
0∈ {(αM ¯q−D)T}+NQad (¯q).
We write for short
Z:= H2(Ω) ×Rl×H2(Ω) and ˜
W:= H2(Ω) ×H2(Ω) ×Rl
and define F:Z→˜
Wand T:Z→2˜
Was follows
F(z) = F(u, q, λ) =
−u+S(−(u·∇)u) + Gc(q)
−λ+S(−(∇¯u)Tλ+ (¯u·∇)λ) + Gw(~el+K(u, p)~ed)
(αMq −D)T
(6.1.1)
for the differentiable part and
T(z) = T(u, q, λ) =
{0}
{0}
NQad (q)
(6.1.2)
for the set-valued part in the generalized equation. Now, we can formulate
the optimality system as
0∈F(¯z) + T(¯z).(6.1.3)
The equation (6.1.3) is a generalized equation. For more details see Alt
[5, 8, 9] Dontchev [33, 35, 36, 37], Goldberg and Tröltzsch [48] and Josephy
[60].
The classical Newton-method is not applicable to this kind of system.
The generalized Newton-method is mentioned in the works above and reads
as
Algorithm 6.1 (GNM).
1. Choose an initial value z0and set k= 0
72
6.1. GENERALIZED EQUATIONS
2. Determine zby the linearized generalized equation
0∈F(zk) + F0(zk)(z−zk) + T(z)(6.1.4)
3. Set k=k+ 1 and zk=z. Goto 2.
Similar to the theory of the classical Newton-method, we need further
assumptions to show the convergence of this generalized Newton-method,
which is closely related to the notion of strong regularity of (6.1.3) and based
on Robinson [82].
Definition 6.2 (Strong regularity).Let ˜z= (˜u, ˜q, ˜
λ)∈Z. The generalized
equation (6.1.3) is said to be strongly regular at ˜z, if there exist open balls
Br1(0) in ˜
Wand Br2(˜z)in Z, with positive constants r1>0,r2>0and
cL>0such that for all perturbations e∈ Br1(0) the linearized equation
e∈F(˜z) + F0(˜z)(z−˜z) + T(z)
admits a unique solution z=z(e)in Br2(˜z)and the Lipschitz property
kz(e1)−z(e2)kZ≤cLke1−e2k˜
W
holds for all e1, e2∈ Br1(0).
In the original paper of Robinson [83], it was assumed that Thas a
closed graph. Dontchev has shown in [35] that this assumption is not needed.
The next theorem provides the possibility to transfer stability results for the
perturbed linearized equation to the perturbed nonlinear equation.
Theorem 6.3. Let ¯z∈Zbe a solution of the generalized equation (6.1.3)
such that this equation is strongly regular at ¯z. Then there exist open balls
Br1(0) and Br2(¯z)such that, for all e∈ Br1(0), the perturbed equation
e∈F(¯z) + T(¯z)
has a unique solution z=z(e)in Br2(¯z)and the solution mapping e7→ z(e)
is Lipschitz-continuous from Br1(0) to Br2(¯z).
Based on the strong regularity, the next theorem states the convergence
of the SQP-method.
Theorem 6.4. Let ¯z= (¯u, ¯p, ¯q, ¯
λ, ¯π)be a solution of (6.1.3) and additionally
let this generalized equation be strongly regular at ¯z. Then there exists an
open ball Br2(¯z)such that for every starting point z1in Br2(¯z)the generalized
Newton method generates a unique sequence {zk}∞
k=1, where zkstays in Br2(¯z)
and we obtain
kzk+1 −¯zkZ≤cGkzk−¯zk2
Z.
with cGindependent of k.
73
CHAPTER 6. CONVERGENCE OF THE SQP-METHOD
For a proof, we refer to [5, 36].
This theorem means that the linearized generalized equation is uniquely
solvable under the suitable assumptions. The sequence of solutions, gen-
erated by the generalized Newton-method, converges quadratically to the
solution of the generalized equation (6.1.3). In the next part we investigate
the solvability of the linearized generalized equation and the relation to the
SQP-method.
Let us consider the linearized equation
0∈F(zi) + F0(zi)(z−zi) + T(z).
This means that we have due to the linearity of K
0∈ −u+S((ui·∇)ui−(ui·∇)u−(u·∇)ui) + Gc(q) + {0}
0∈ −λ+S(−(∇u)Tλi+ (u·∇)λi−(∇ui)Tλ+ (ui·∇)λ
+ (∇ui)Tλi−(ui·∇)λi) + Gw(el+˜
K(u, p)ed) + {0}
0∈(αMq −D(zi)−D0(zi)(z−zi))T+NQad (q).
(6.1.5)
The first relation of (6.1.5) is equivalent to the boundary value problem
−ν∆u+ (ui·∇)u+ (u·∇)ui+∇p= (ui·∇)uiin Ω
div u= 0 in Ω
u=
l
X
i=1
eiqion Γc,
u= 0 on Γ\Γc
(6.1.6)
and the second to the problem
−ν∆λ+∇π=−(∇u)Tλi−(∇ui)T(λ−λi)
+ (ui·∇)(λ−λi)+(u·∇)λiin Ω
div λ= 0 in Ω
λ=~el+˜
K(u, p)~edon Γw
λ= 0 on Γ\Γwall.
(6.1.7)
Considering NQad (q), we additionally get
(αMq −D(zi)−D0(zi)(z−zi))T(˜q−q)≥0∀˜q∈Qad.(6.1.8)
We see that the linearized generalized equation (6.1.4) is equivalent to (6.1.6)-
(6.1.8) which is similar to the first-order necessary optimality system, see
74
6.1. GENERALIZED EQUATIONS
(6.1.9). In reality, this is the necessary optimality system for the linear-
quadratic problem of the SQP-method, mentioned above. The cost functional
of the linear-quadratic problem was defined by
Jk(u, q) = ∇J(uk, qk)(u−uk, q −qk) + 1
2Lzz(uk, qk, λk)[(u−uk, q −qk)]2.
After expanding, it is with g=Hc(q) = Pl
i=1 eiqiequal to
Jk(u, q) = Z
Γwν∂(u−uk)
∂n−(p−pk)n·~elds
+K(uk, pk)
Z
Γwν∂(u−uk)
∂n−(p−pk)n·~edds
+αZ
Γc
gk(g−gk)ds+b(u−uk, u −uk, λk)
+˜
K(u−uk, p −pk)
Z
Γwν∂(u−uk)
∂n−(p−pk)n·~edds
=Z
Γwν∂(u−uk)
∂n−(p−pk)n·~elds
+K(u, p)
Z
Γwν∂(u−uk)
∂n−(p−pk)n·~edds
+αZ
Γc
gk(g−gk)ds+b(u−uk, u −uk, λk).
Now, we have shown the connection between the generalized Newton-method
and the SQP-method. The SQP-method is a kind of a generalized Newton-
method to solve the optimality system
−ν∆u+ (u·∇)u+∇p=fin Ω
div u= 0 in Ω
u=
l
X
i=1
eiqion Γc,
u= 0 on Γ\Γc.
(6.1.9a)
75
CHAPTER 6. CONVERGENCE OF THE SQP-METHOD
−ν∆λ+ (∇¯u)Tλ−(¯u·∇)λ+∇π= 0 in Ω
div λ= 0 in Ω
λ=~el+K(¯u, ¯p)~edon Γw
λ= 0 on Γ\Γw,
(6.1.9b)
(αM ¯q−D)T(q−¯q)≥0∀q∈Qad (6.1.9c)
obtained analogously to (4.3.4a)-(4.3.4c).
The optimality system of the linear-quadratic system is equivalent to the
linearized generalized system (6.1.4). For the finite-dimensional case, one can
find an example in [42] and for the case of infinite-dimensional optimization,
see [6] and [33]. In the next sections, we want to prove the convergence of
the SQP-method based on Theorem 6.4. That means that we have to show
the assumptions in this theorem.
6.2 Perturbed optimization problem
There is a lot of literature about convergence of the SQP-method. Alt [6, 7],
Alt and Malanowski [10, 11], Malanowski [68, 67] and Dontchev et. al. [34]
proved convergence for optimal control problems subject to ODEs. Optimal
control problems related to PDEs were investigated in Hinze and Hinter-
müller [53], Goldberg and Tröltzsch [48], Kupfer and Sachs [64], Tröltzsch
[100, 101], and Volkwein [107].
In the last section, we have shown that the SQP-method for the optimal
boundary control problem is a variation of the generalized Newton-method
for solving the first-order necessary optimality conditions. In this section of
this thesis, we follow the work of Unger [104] and Wachsmuth [109].
In the next step, we want to verify continuous differentiability of F, Lip-
schitz continuity of F0, and strong regularity of the optimality system to use
Theorem 6.4 and to show convergence of the SQP-method.
Theorem 6.5. The function Fdefined by (6.1.1) is continuous differentiable
and the derivative B0(u)is Lipschitz continuous.
The proof is analogous to [109], Corollary 5.2.
Now, let us come to the more complex assumption of strong regular-
ity. Therefore, we want to consider the pertubed linearized at ¯zgeneralized
equation
e∈F(¯z) + F0(¯z)(z−¯z) + T(z)(6.2.1)
as an optimality system of another optimal control problem. Let
¯z= (¯u, ¯p, ¯q, ¯
λ, ¯π)fulfill the the first-order necessary optimality conditions of
76
6.2. PERTURBED OPTIMIZATION PROBLEM
(PI)and because of this be a solution of the generalized equation (6.1.3).
The perturbed equation (6.2.1) with e= (˜eu,˜eλ,˜eq)∈ˆ
Wat ¯zis equivalent
to
0∈ −˜eu−u+S(−(¯u·∇)u−(u·∇)¯u+ (¯u·∇)¯u) + Gc(q) + {0}
0∈ −˜eλ−λ+S(−(∇u)T¯
λ+ (u·∇)¯
λ−(∇¯u)Tλ+ (¯u·∇)λ
+ (∇¯u)T¯
λ−(¯u·∇)¯
λ) + Gw(el+˜
K(u, p)ed) + {0}
0∈(−˜eq+αMq −D(¯z)−D0(¯z)(z−¯z))T+NQad (q),
with ˆu=u+ ˜eu,ˆp=p, ˆq=q, ˆ
λ=λ+ ˜eλand ˆπ=πequivalent to
0∈ −ˆu+S(−(¯u·∇)(ˆu−˜eu)−((ˆu−˜eu)·∇)¯u+ (¯u·∇)¯u) + Gc(ˆq) + {0}
0∈ −ˆ
λ+S(−(∇(ˆu−˜eu))T¯
λ+ ((ˆu−˜eu)·∇)¯
λ+ (∇¯u)T¯
λ−(¯u·∇)¯
λ
−(∇¯u)T(ˆ
λ−˜eλ) + (¯u·∇)(ˆ
λ−˜eλ))
+Gw(~el+˜
K(ˆu−˜eu, p)~ed) + {0}
0∈(−˜eq+αM ˆq−D(¯z)−D0(¯z)((ˆu−˜eu,ˆp, ˆq, ˆ
λ−˜eλ,ˆπ)−¯z))T+NQad (q).
The last system is with
eu= (¯u·∇)˜eu+ (˜eu·∇)¯u,
e1
λ= (∇˜eu)T¯
λ−(˜eu·∇)¯
λ+ (∇¯u)T˜eλ−(¯u·∇)˜eλ,
e2
λ=−K(˜eu,0),
Neq=−˜eq−˜
D
(6.2.2)
and
˜
Di=ZΓc
(ν∂˜eλ
∂n+ (¯u·n)˜eλ+ (˜eu·n)¯
λ+ (¯u·n)¯
λ)eids,
ˆ
Di=ZΓc
(ν∂ˆ
λ
∂n−ˆπn+ (¯u·n)ˆ
λ+ (ˆu·n)¯
λ)eids
N∈Nl, Ni=ZΓc
eids
for i= 1, . . . , l equivalent to
0∈ −ˆu+S(−(¯u·∇)ˆu−(ˆu·∇)¯u+ (¯u·∇)¯u) + Gc(ˆq) + S(eu) + {0}
0∈ −ˆ
λ+S(e1
λ) + Gw(e2
λ) + Gw(~el+˜
K(ˆu, p)~ed) + {0}
+S(−(∇ˆu)T¯
λ+ (ˆu·∇)¯
λ−(∇¯u)Tˆ
λ+ (¯u·∇)ˆ
λ+ (∇¯u)T¯
λ−(¯u·∇)¯
λ)
0∈(Neq+αM ˆq−ˆ
D)T+NQad (q).
(6.2.3)
77
CHAPTER 6. CONVERGENCE OF THE SQP-METHOD
The system (6.2.3) corresponds to the following system of the perturbed
state equation
−ν∆ˆu+ (¯u·∇)ˆu+ (ˆu·∇)¯u+∇p= (¯u·∇)¯u+euin Ω
div ˆu= 0 in Ω
ˆu=
l
X
i=1
eiˆqion Γc,
ˆu= 0 on Γ\Γc,
(6.2.4)
the perturbed adjoint equation
−ν∆ˆ
λ+ (∇¯u)Tˆ
λ−(¯u·∇)ˆ
λ+∇π=−(∇ˆu)T¯
λ+ (ˆu·∇)¯
λ+e1
λ
+ (∇¯u)T¯
λ−(¯u·∇)¯
λin Ω
div λ= 0 in Ω
λ=~el+˜
K(ˆu, p)~ed+e2
λon Γw
λ= 0 on Γ\Γw
(6.2.5)
and the associated variational inequality
(αM ˆq−ˆ
D+Neq)T(q−¯q)≥0∀q∈Qad.(6.2.6)
Now, we see that this is the optimality problem of the following pertubated
linear-quadratic optimization problem (Pe), where we write (u, p, q, λ, π) =
(ˆu, ˆp, ˆq, ˆ
λ, ˆπ)for simplicity:
min Je(u, q) := Z
Γw∂u
∂n−pn·~elds+c
2
Z
Γw∂u
∂n−pn·~edds−D0
2
+Z
Γc
|Hc(q)|2ds+Z
Ω
euu+e1
λλdx+Z
Γc
eqHc(q)ds
+Z
Γw
e2
λλds+b(u−¯u, u −¯u, ¯
λ)
subject to the pertubated Navier-Stokes equations
−ν∆u+ (¯u·∇)u+ (u·∇)¯u+∇p= (¯u·∇)¯u+euin Ω
div u= 0 in Ω
u=
l
X
i=1
eiqion Γc,
u= 0 on Γ\Γc,
(6.2.7)
78
6.3. A MODIFIED PROBLEM
and the control constraints
q∈Qad.
Let us now consider e:= (eu, e1
λ, e2
λ, eq)∈W:= H2(Ω) ×H2(Ω) ×R×Rl.
Then, we obtain due to (6.2.2)
kekW≤Ck˜ek˜
W
with kekW=keukH2(Ω) +ke1
λkH2(Ω) +|e2
λ|+|eq|,
kekW=k˜eukH2(Ω) +k˜eλkH2(Ω) +|˜eq|.
6.3 A modified problem
Because there is no convexity of the cost functional Jein directions not
included in (SSC), we need
˜
Qad := {˜q∈Qad : ˜q= ¯qon A}
to not allow changes of the control function qon the active sets. Let us
denote the new problem by (˜
Pe)consisting of the cost functional Je, the
state equation (6.2.7) and the control constraint q∈˜
Qad. Finding a solution
for the problem (˜
Pe)is equivalent to finding a solution for the linearized and
perturbed generalized equation (6.2.1), with
T(z) = ({0},{0},N˜
Qad (q))T.
First, we investigate the existence of a solution for the perturbed linear-
quadratic optimal control problem (˜
Pe). Then, we want to prove that the
optimal control of (˜
Pe)is Lipschitz-continuous with respect to the perturba-
tion e. The idea is to investigate strong regularity of the generalized equation
e∈F(z) + (0,0,N˜
Qad (q)) (6.3.1)
and then transfer the result to the original problem.
6.3.1 Existence of a solution
The next theorem proves the solvability of (˜
Pe).
Theorem 6.6. Assume that ¯z= (¯u, ¯p, ¯q, ¯
λ, ¯π)satisfy the optimality system
and the coercivity condition (SSC00)at the begin of this chapter. Then (˜
Pe)
admits a unique optimal control ¯ue.
79
CHAPTER 6. CONVERGENCE OF THE SQP-METHOD
Proof. Let c∈Rbe a generic constant. Denoting the Lagrangian belonging
to (˜
Pe)by Le, we derive
Le
qq(q) = Lqq(¯q)(6.3.2)
for all q, because the perturbation appears only linear.
Taking q1, q2∈˜
Qad with associating u1, u2of (6.2.7), the pair (u1−u2, q1−
q2)fits to the assumption of (SSC00)and we obtain
(q1−q2)Le
qq(q)(q1−q2)=(q1−q2)Lqq(¯q)(q1−q2)>0.
for q16=q2.Because of this, the problem (˜
Pe)is convex on ˜
Qad. Thus, (˜
Pe)
is uniquely solvable as a linear-quadratic optimization problem with strongly
convex cost functional with modifications described in Section 3.3. We denote
the unique solution by ¯ue.
For more details, see [104].
6.3.2 Lipschitz stability
We still need to prove the Lipschitz-continuity of the perturbation to solution
mapping e7→ (ue, ge, λe)of (˜
Pe)to show strong regularity.
Theorem 6.7. Let (¯z)fulfill the coercivity condition (SSC00). Then the
solution mapping e:= (eu, e1
λ, e2
λ, eq)7→ ze= (ue, pe, qe, λe, πe)is Lipschitz-
continuous from Wto Z.
Proof. Let z1, z2be two elements of Zand qi,i= 1,2be the optimal control
functions of the optimization problem (˜
Pe)with the associated states uiand
adjoints λi. Let us furthermore define the differences z:= z1−z2,q:= q1−q2,
u:= u1−u2and λ:= λ1−λ2.
So, the variational inequality with the constraint qi∈˜
Qad looks as
(αMqi−ˆ
D(ui, λi, πi)−Neq,i)T(q−qi)≥0,∀q∈˜
Qad.(6.3.3)
Testing this inequality with q2for i= 1 and q1for i= 2 and adding both,
we derive ˆ
D(u, λ, π)Tq+ (Neq)Tq≥(αMq)Tq. (6.3.4)
The difference uis the weak solution of the state equation
−ν∆u+ (¯u·∇)u+ (u·∇)¯u+∇p=euin Ω
div u= 0 in Ω
u=g:=
l
X
j=1
ejqjon Γc,
u= 0 on Γ\Γc,
(6.3.5)
80
6.3. A MODIFIED PROBLEM
Testing (6.3.5) with (λ, π)=(λ1−λ2, π1−π2), we obtain
a(u, λ)−(ν∂u
∂n−pn, λ)L2(Γ) + (p, div λ)L2(Ω)
+h(¯u·∇)u+ (u·∇)¯u, λi(H1(Ω))0,H1(Ω) = (eu, λ)L2(Ω)
(div u, π)L2(Ω) = 0
τΓcu=g
(6.3.6)
and testing the adjoint
−ν∆λ+ (∇¯u)Tλ−(¯u·∇)λ+∇π=−(∇u)T¯
λ+ (u·∇)¯
λ+e1
λin Ω
div λ= 0 in Ω
λ=~el+˜
K(u, p)~ed+e2
λon Γw
λ= 0 on Γ\Γw
(6.3.7)
with (u, p) = (u1−u2, p1−p2), we get
a(λ, u)−(ν∂λ
∂n−πn+ (¯u·n)λ+ (u·n)¯
λ, u)L2(Γ)
+h(∇¯u)Tλ−(¯u·∇)λ, ui(H1(Ω))0,H1(Ω)
+h(∇u)T¯
λ−(u·∇)¯
λ, ui(H1(Ω))0,H1(Ω) + (π, div u)L2(Ω) = (e1
λ, u)L2(Ω)
(6.3.8a)
(div λ, p)L2(Ω) = 0 (6.3.8b)
τΓwλ=~el+˜
K(u, p)~ed+e2
λ.(6.3.8c)
The systems (6.3.6) and (6.3.8) are equivalent to
a(u, λ)−ν(∂u
∂n−pn,~el+˜
K(u, p)~ed+e2
λ)L2(Γw)
+h(¯u·∇)u+ (u·∇)¯u, λi(H1(Ω))0,H1(Ω) = (eu, λ)L2(Ω)
and
a(λ, u)−ν(∂λ
∂n−πn+ (¯u·n)λ+ (u·n)¯
λ, g)L2(Γc)
+h(¯u·∇)u+ (u·∇)¯u, λi(H1(Ω))0,H1(Ω)
+h(∇u)T¯
λ−(u·∇)¯
λ, ui(H1(Ω))0,H1(Ω) = (e1
λ, u)L2(Ω).
Considering them with
h(∇u)T¯
λ−(u·∇)¯
λ, ui(H1(Ω))0,H1(Ω) =h(u·∇)u, ¯
λi(H1(Ω))0,H1(Ω)
81
CHAPTER 6. CONVERGENCE OF THE SQP-METHOD
we obtain the equality
(eu, λ)L2(Ω)+(ν∂u
∂n−pn,~el+˜
K(u, p)~ed+e2
λ)L2(Γw)
= (e1
λ, u)L2(Ω) −h(u·∇)u, ¯
λi(H1(Ω))0,H1(Ω)
+ (ν∂λ
∂n−πn+ (¯u·n)λ+ (u·n)¯
λ, g)L2(Γc)
or equivalent
h(u·∇)u, ¯
λi(H1(Ω))0,H1(Ω) −˜
K(u, p)(ν∂u
∂n−pn,~ed)L2(Γw)
= (e1
λ, u)L2(Ω) −(eu, λ)L2(Ω) −(ν∂u
∂n−pn, e2
λ)L2(Γw)
+ (ν∂λ
∂n−πn+ (¯u·n)λ+ (u·n)¯
λ, g)L2(Γc)
−(ν∂u
∂n−pn,~el)L2(Γw).
(6.3.9)
Let us additionally consider ˜uas a weak solution of (6.3.5) with q= 0. Thus,
(u−˜u, q)fits to the assumption of the coercivity condition (SSC00), due to
qi= (qe)1,i −(qe)2,i = (¯qe)i−(¯qe)i= 0 on Aby (qe)1,(qe)2∈˜
Qad and (u−˜u)
is the solution of the associated at ¯qlinearized state equation.
Furthermore, we obtain
0<Lvv(¯z)[u−˜u, q]2
=Lqq(¯z)[q]2
| {z }
∗1
+Luu(¯z)[u]2
| {z }
∗2
−2Luu(¯z)[u, ˜u]2
| {z }
∗3
+Luu(¯z)[˜u]2
| {z }
∗4
.(6.3.10)
Adding the terms ∗1and ∗2leads to
Lqq(¯z)[q]2+Luu(¯z)[u]2=αkgk2
L2(Γc)−h(u·∇)u, ¯
λi(H1(Ω))0,H1(Ω)
+˜
K(u, p)(ν∂u
∂n−pn,~ed)L2(Γw)
and with (6.3.9) and kgkL2(Γc)= (Mq)Tqto
Lqq(¯z)[q]2+Luu(¯z)[u]2=α(Mq)Tq+ (eu, λ)L2(Ω) −(e1
λ, u)L2(Ω)
+ (ν∂u
∂n−pn, e2
λ+~el)L2(Γw)
−(ν∂λ
∂n−πn+ (¯u·n)λ+ (u·n)¯
λ, g)L2(Γc).
82
6.3. A MODIFIED PROBLEM
The inequality (6.3.4) leads to
Lqq(¯z)[q]2+Luu(¯z)[u]2≤(eu, λ)L2(Ω) −(e1
λ, u)L2(Ω) + (Neg)Tq
+ (ν∂u
∂n−pn, e2
λ+~el)L2(Γw)+ˆ
D(u, λ, π)Tq
−(ν∂λ
∂n−πn+ (¯u·n)λ+ (u·n)¯
λ, g)L2(Γc)
≤ ke1
λkL2(Ω)kukL2(Ω) +keukL2(Ω)kλkL2(Ω) + (Neg)Tq
+kν∂u
∂n−pnkL2(Γc)(ke2
λkL2(Γc)+k~elkL2(Γc))
−(ν∂λ
∂n−πn+ (¯u·n)λ+ (u·n)¯
λ, g)L2(Γc)
+ (ν∂λ
∂n−πn+ (¯u·n)λ+ (u·n)¯
λ, g)L2(Γc)
≤ kekW(kukH2(Ω) +kpkH1(Ω) +kλkH2(Ω))+(Neg)Tq.
(6.3.11)
We get
k˜ukH2(Ω) +k˜pkH1(Ω) ≤ckeukH2(Ω) ≤ckekW
because ˜uis the weak solution of the at (¯u, ¯g)linearized Navier-Stokes equa-
tions, similar to the regularity Assumption of non-singularity 2.25. Another
assertion is that we obtain by Theorem 2.11 that if (u, p)is the solution of the
nonhomogenous Stokes system (2.2.1) with a sufficiently smooth boundary,
then we get
kukH2(Ω) +kpkH1(Ω) ≤c(kfkL2(Ω) +kgkH3/2(Γc)).
Transfering this over to the linearized Navier-Stokes equations (6.3.5) with
g= 0, we obtain
kukH2(Ω) +kpkH1(Ω) ≤c(k−(¯u·∇)u−(u·∇)¯u+eukL2(Ω))
≤c(k(¯u·∇)u+ (u·∇)¯ukL2(Ω) +keukL2(Ω))
≤c(kukH2(Ω) +keukL2(Ω))
≤c(keukL2(Ω))
and this also leads to
k˜ukH2(Ω) +k˜pkH1(Ω) ≤ckekW.
83
CHAPTER 6. CONVERGENCE OF THE SQP-METHOD
Adding the terms ∗3and ∗4of (6.3.10), we obtain
|2Luu(¯z)[u, ˜u]2|+|Luu(¯z)[˜u]2| ≤ 2|h(u·∇)˜u+ (˜u·∇)u, ¯
λi(H1(Ω))0,H1(Ω)|
+|h2(˜u·∇)˜u, ¯
λi(H1(Ω))0,H1(Ω)|
+˜
K(˜u, ˜p)( ˜
K(˜u, ˜p) + ˜
K(u, p))
≤2kuk1/2
L2(Ω)kuk1/2
H2(Ω)k˜ukH1(Ω)k¯
λkL2(Ω)
+ck˜ukH2(Ω)k~edkL2(Ω) +ck˜pkH2(Ω)k~edkL2(Ω)
+k˜uk1/2
L2(Ω)k˜uk1/2
H2(Ω)k˜ukH1(Ω)k¯
λkL2(Ω)
+ (kukH2(Ω) +kpkH1(Ω))(k˜ukH2(Ω) +k˜pkH1(Ω))
+ (k˜ukH2(Ω) +k˜pkH1(Ω))(k˜ukH2(Ω) +k˜pkH1(Ω))
(6.3.12)
with similar arguments to Lemma 2.14. Furthermore, this leads to
|2Luu(¯z)[u, ˜u]2|+|Luu(¯z)[˜u]2| ≤ c((kukH2(Ω) +kpkH1(Ω) +kekW)kekW+kekW)
≤ckekW(kukH2(Ω) +kpkH1(Ω) +kekW+ 1).
(6.3.13)
Let us now summarize (6.3.10)-(6.3.13) to:
0< ckekW(kekW+kukH2(Ω) +kpkH1(Ω)
+kλkH2(Ω) +kπkH1(Ω))+(Neq)Tq. (6.3.14)
Because (6.3.5) and (6.3.7) are at (¯u, ¯g)linearized Navier-Stokes equations,
we obtain from the non-singularity assumption, similar to the Definition 2.25:
kukH2(Ω) +kpkH1(Ω) ≤c(kgkL2(Γc)+kekW)≤c(|q|+kekW),
kλkH2(Ω) +kπkH1(Ω) ≤c(k~el+˜
K(u, p)~ed+e2
λkL2(Γw)+kukH2(Ω) +ke1
λkL2(Ω))
≤c(kukH2(Ω) +kpkH1(Ω) +ke1
λkL2(Ω) + +ke2
λkL2(Γw))
≤c(kgkL2(Γc)+kekW)≤c(|q|+kekW).
(6.3.15)
The inequality (6.3.14) yields to
0< ckekW(kekW+kgkL2(Γc))
≤ckek2
W+δ
2kgk2
L2(Γc)≤c(|q|+kekW)
and proves Lipschitz-continuity of e7→ qfrom Wto Rl. All this leads to a
Lipschitz continuity of e7→ (u, q, λ)from Wto Z=H2(Ω)×Rl×H2(Ω).
84
6.4. STRONG REGULARITY OF THE ORIGINAL PERTURBED
PROBLEM
6.4 Strong regularity of the original perturbed
problem
Now, we return to the original perturbed problem (Pe). To prove strong reg-
ularity of the generalized equation (6.1.3), we have to search for a solution
qe=q(e)in Qad. The first idea is to take qe∈˜
Qad, solving (˜
Pe). We will
show in this section that qeis also a solution of (Pe)for a sufficiently small
perturbation e∈W. Therefore, we have to investigate the optimal control
qeon the active set A. In this section, we closely follow again Wachsmuth
and want to refer to [109, Chapter 5, Section 4], for the ideas of the proofs
of the following theorems.
Lemma 6.8. Let ¯z= (¯u, ¯p, ¯q, ¯
λ, ¯π)fulfill the coercivity condition (SSC00)and
let qa, qb∈Rl. Then exists ρe>0and σ > 0such that the optimal control
function qeof (˜
Pe)is active for all e∈Wwith kekW≤ρe. This means
(αMqe−ˆ
D+Neq)T
i>σ
2on A+,
(αMqe−ˆ
D+Neq)T
i<−σ
2on A−
(6.4.1)
and sign(αMq −ˆ
D+Neq)T
i=sign(αM ¯q−ˆ
D)T
ion A.
Proof. Theorem 6.7 garantuees that the mapping e7→ (ue, pe, qe, λe, πe)is
Lipschitz continuous from Wto Z. Furthermore, it is easy to see that ˆ
D(z) =
RΓc(ν∂λ/∂n−πn+ (¯u·n)λ+ (u·n)¯
λ)eidsis Lipschitz continuous from
L2(Ω) ×L2(Ω) ×Rl×L2(Ω) ×L2(Ω) to Rl. Due to qe∈Rl,there exists a
σ > 0with
|αM ¯q−ˆ
D(z)|> σ.
Considering i∈ A+,we obtain
σ < (αM ¯q−ˆ
D(z))T
i
= (αM ¯q−ˆ
D(z))T
i−(αM ¯q−ˆ
D(z) + Neq)T
i+ (αM ¯q−ˆ
D(z) + Neq)T
i
≤ckekW+ (αM ¯q−ˆ
D(z) + Neq)T
i.
Taking ρesufficiently small yields
(αM ¯q−ˆ
D(z) + Neq)T
i>σ
2
Analogously leads i∈ A−to
(αM ¯q−ˆ
D(z) + Neq)T
i<−σ
2.
85
CHAPTER 6. CONVERGENCE OF THE SQP-METHOD
An important consequence is that the control function qeeven satisfies
the variational inequality based on the admissible set Qad and not only for
˜
Qad ⊂Qad.
Lemma 6.9. With the assumptions and ρeas in the last Lemma 6.8, the
control function qeassociated to a perturbation e, fulfilling kekW≤ρe, satis-
fies
(αMqe−ˆ
D+Neq)T(q−qe)≥0∀q∈Qad.
Proof. Let q∈Qad.Then we are able to split (αMqe−ˆ
D+Neq)T(q−qe)
into
(αMqe−ˆ
D+Neq)T(q−qe) = X
i∈A+
(αMqe−ˆ
D+Neq)T
i(q−qe)i
+X
i∈A−
(αMqe−ˆ
D+Neq)T
i(q−qe)i
+X
i/∈A
(αMqe−ˆ
D+Neq)T
i(q−qe)i
with qe∈˜
Qad.The third sum is nonnegative due to the fact that it is part
of the optimality system of (˜
Pe).
Additionally, we have
qe= ¯q=qaon A+,
qe= ¯q=qbon A−.
Because of kekW≤ρe,Lemma 6.8 leads to
(αMqe−ˆ
D+Neq)T
i>σ
2on A+,
(αMqe−ˆ
D+Neq)T
i<−σ
2on A−.
Thus, we obtain
(αMqe−ˆ
D+Neq)T(q−qe)≥0∀q∈Qad.
Remark 6.10. The triple (ue, ge, λe)of Theorem 6.9 fulfills the optimality
system of the optimization perturbed problem (Pe)which is equivalent to the
linearized and perturbed generalized equation (6.2.1).
The next theorem shows with the help of the equality (6.3.2) that (ue, qe, λe)
minimizes the optimization problem (Pe).
86
6.4. STRONG REGULARITY OF THE ORIGINAL PERTURBED
PROBLEM
Theorem 6.11. Under the same assumptions as in the two last lemma, there
exist ρe, ρq>0such that the control qebelonging to a perturbation e∈W
satisfying kekW≤ρeis locally optimal for the optimal control problem (Pe)
and fulfills
Je(ue, qe)≤Je(u, q)
for all q∈Qad fulfilling |q−qe| ≤ ρq, where uand ueare the weak solutions
of (6.2.4) associated to qand qe, respectively.
The proof is very similar to [109, Theorem 5.15].
Now, we have shown the existence of a solution of the linearized and
perturbed equation (6.2.1).
Theorem 6.11 shows that qeis the unique optimal solution of Problem
(Pe)in Bρq(¯q)with perturbations ein Bρe(0). By Theorem 6.7 uqand λq
are in Bcuρe(¯u)and Bcλρe(¯
λ)with the Lipschitz-constants cuand cλgiven by
Theorem 6.7. This leads to the unique solvability of (6.2.1) in Bcuρe(¯u)×
Bρq(¯q)×Bcλρe(¯
λ)for perturbations ein Bρe(0).
This yields the strong regularity of the generalized equation (6.1.3).
The investigations of this chapter have shown that we only find a local
solution of the linearized subproblems of the SQP-method in a neighborhood
of the reference solution ¯v= (¯u, ¯q). The idea is now to modify Qad to
Qρ
ad := Qad ∩{h∈Rl:|h−¯q| ≤ ρ}
to have the solution of the linearized subproblems in a close neighborhood of
the reference solution. See [104] for more details.
Altogether, it follows the local convergence of the SQP-method, see The-
orem 6.4.
Theorem 6.12. Let ¯z∈Zfulfill the coercivity condition (SSC00). Then
exist ρ > 0such that the SQP-method with control constraint Qρ
ad generates
a uniquely determined sequence (uk, qk, λk),qk∈Qρ
ad for every starting point
(u0, q0, λ0), with q0∈Qρ
ad and we obtain
kuk+1 −¯ukH2(Ω) +|qk+1 −¯q|+kλk+1 −¯
λkH2(Ω)
≤c(kuk+1 −¯uk2
H2(Ω) +|qk+1 −¯q|2+kλk+1 −¯
λk2
H2(Ω)).
87
CHAPTER 6. CONVERGENCE OF THE SQP-METHOD
88
Chapter 7
The nonstationary case
After investigating the steady-state problem, we want to focus in this chap-
ter on the optimization problem subject to the nonstationary Navier-Stokes
equations. A simplified model, similar to the stationary case, could look as
follows. We want to minimize the negative lift
min
u,g J(u, g) := −ZZ
Σν∂u
∂n−pn·~eldxdt+α
2kgkL2(Σ)
subject to the nonstationary Navier-Stokes equations describing the motion
of the fluid around the body
ut+ν∆u+ (u·∇)u+∇p=fon Q,
div u= 0 on Q,
u=gon Σ,
u(0) = u0on Ω
and the control constraint
g∈Gad,
where Q:= Ω ×(0, T)with its boundary Σ := ∂Q = Γ ×(0, T)with a
fixed time T. But in contrast to the stationary case, we allow high Reynolds
numbers in the nonstationary situation hand consider a problem closer to
the real setting of the high-lift configuration, see e.g. [19, 89, 98].
In this situation, we have to deal with turbulence’s, which we simulated
by a k-ω-Wilcox98 model, see [113], where the equations for kand ωare
89
CHAPTER 7. THE NONSTATIONARY CASE
given by:
u∂u
∂x +v∂u
∂y =∂
∂y[(ν+νT)∂u
∂y ]
u∂k
∂x +v∂k
∂y =νT(∂u
∂y )2−β∗ωk +∂
∂y[(ν+σ∗νT)∂k
∂y ]
u∂ω
∂x +v∂ω
∂y =αω
kνT(∂u
∂y )2−βω2+∂
∂y[(ν+σ∗νT)∂ω
∂y ]
νT=α∗k
ω,
(7.0.2)
where uand vare velocity components in the streamwise xand normal y
directions, νis the kinematic molecular viscosity, νTis the kinematic eddy
viscosity and β, β∗, σ, σ∗are parameters, which are defined in [113].
Due to the high dimension of the discretized equations, the computing
times for any forward solution of the model are extremely large so that a
mathematical optimization of the periodic actuation is fairly unrealistic. In
[19], a generic high-lift configuration was investigated and one forward so-
lution took about 48 hours. In the case of the SCCH configuration, the
computation time was nearly twice that number.
Let us mention that in this nonstationary case, we want to consider a
special example with the following setting.
Setting 7.1. We consider the incompressible two-dimensional flow over the
swept constant chord half (SCCH) high-lift configuration, see Figure 7.1. The
chord length cis denoted by Lref = 1.275 and the inflow by u∞= 1. The
chord length is the length of the wing in the flow direction. The Reynolds
number is Re =u∞c/ν, where νis the kinematic viscosity of the fluid. The
leading edge slat deflection angle is 26.5°, the flap deflection is 37°and the
angle of attack of the wing is 6°. The periodic actuation is introduced by a
zero-net-mass-flux actuator on a small slit on the flap, where the flow fully
separates. The actuation velocity is
g(t) = Bcos(Ωat),(7.0.3)
where Ωa= 2πStais the angular actuation frequency, Bthe actuation ampli-
tude, Sta=fac/u∞the Strouhal number and fais the actuation frequency.
Analogously, we define Stn=fnc/u∞with the vortex-shedding frequency fn.
The actuation intensity is characterized by the dimensionless coefficient
Cµ=H
cB
u∞2
90
7.1. MODEL REDUCTION
with the slot width H= 0.001238cfl and the relative chord length cfl = 0.254c.
The full k−ωWilcox98 model was solved by unsteady Reynolds-averaged
Navier-Stokes (URANS) equations with the ELAN code1. The with cfl non-
dimensionalized natural Strouhal number is Stn
fl =fncfl/u∞= 0.32. The
actuation is described by a momentum coefficient of Cµ= 405 ×10−5and an
actuation frequency of Sta
fl = 0.6.
Figure 7.1: The SCCH high-lift configuration, where the periodic excitation
is implemented on the flap.
We think that a model reduction is advisable, because of the reasons
above. Our goal is to establish a reduced-order model (ROM) as a basis for
our optimization problem.
7.1 Model reduction
The topic of model reduction is currently in great demand by engineers. A
widely used method is POD [4, 56, 62, 63, 111]. In the case of the high-lift
configuration, the application of standard POD does not align to the target
of robust dynamical least-order models for the real flow. To establish the
1Developed at the Computational Fluid Dynamics and Aeroacoustics Group (Professor
F. Thiele) at the TU Berlin.
91
CHAPTER 7. THE NONSTATIONARY CASE
reduced-order model (ROM), the computed POD basis has to be inserted as
a Galerkin basis in the full Wilcox98 model, see (7.0.2). The associated
implementation would be a time consuming task.
There were several approaches to deal with these problems, e.g. an ex-
tension of POD to data compression of multiple operation points, see [59]
for sequential POD or [93] for DPOD. We follow an alternative approach
suggested in [66, 74] of a canonical reduction with parameter identification.
Here, a very small system of nonlinear ODEs is adapted to the previously
computed flows in the actuated and nonactuated case. This small system is
easily tractable by optimization.
In the next chapter, we will go into details of this technique, present
our modification and report about first experiences in the simplified two-
dimensional Setting 7.1. Our numerical results are promising for future op-
timization tasks.
We also refer to [96] for a similar approach.
Let us first briefly introduce the proper orthogonal decomposition (POD),
which we will need to establish the reduced-order model.
7.2 Proper orthogonal decomposition POD
This section is based on the theory in [105, 106, 108, 107, 111]. In this section,
we want to introduce the proper orthogonal decomposition (POD) and the
way to calculate the POD basis by minimizing a least-square error formula.
There are two cases, the infinite-dimensional and the finite-dimensional
one, of the POD basis. We will consider the finite-dimensional one, because
we want to concentrate on real computations and there we don’t have the
whole trajectory u(t). Therefore, we get an ensemble of snapshots. After-
wards, it is possible to prove that the POD basis is the best orthogonal
system in the ensemble capturing more kinetic energy than any other one
having the same basis number.
So, let ˆu(ti)∈Vbe the Nsnapshots computed by the full dynamical
system with the Setting 7.1 at given times tiwith i= 1, . . . , N,0 = t1≤
t2≤. . . ≤tN−1≤tN=Tand N∈Nand at least one snapshot has to be
non-zero. In the next chapter, we apply the POD method for both cases, the
natural and the actuated one. Let us just explain the method for one case.
We define ˆui:= ˆu(ti),i= 1, . . . , N.
Furthermore, we define VNas the span of the Nsnapshots
VN:= span{ˆu(t1),...,ˆu(tN)}
92
7.2. PROPER ORTHOGONAL DECOMPOSITION POD
with 1≤dimVN≤N.
Let {Φk}N
k=1 denotes an orthonormal basis for VN, which has still to be
computed, then each snapshot ˆu(ti),i= 1, . . . , N, can be expressed by
ˆu(ti) =
N
X
k=1hˆu(ti),ΦkiVΦkfor i= 1, . . . , N.
The idea is to expect that only some of the orthonormal basis func-
tions {Φk}N
k=1 keep most of the kinetic energy so that they can represent
the structure of the snapshots as good as possible. Let therefore M∈Rwith
0< M ≤Nbe given.
Mathematically, we can formulate the problem of finding the orthonormal
system {Φk}M
k=1 by
min
Φk
N
X
j=1
αjkˆu(tj)−
M
X
k=1hˆu(tj),ΦkiΦkk2
V
subject to hΦi,ΦjiV=δij,1≤i, j ≤M
(7.2.1)
where αi’s stand for the trapezoidal weights
α1=t2−t1
2, αi=ti+1 −ti−1
2for 2≤i≤N−1, αn=tN−tN−1
2.
For the next remark, see [111] Chapter 3.
Remark 7.2.
•The trapezoidal approximation for the integral
I(u) =
T
Z0
kˆu(t)−
M
X
k=1hˆu(t),ΦkiΦkk2
Vdt
is
In(u) =
N
X
j=1
αjkˆu(tj)−
M
X
k=1hˆu(tj),ΦkiΦkk2
V
for all u∈C([0, T], V )and it follows that lim
n→∞ In(u) = I(u).
•The least-square problem is equivalent to the largest mean square pro-
jection of the snapshot, namely
max
Φk
N
X
j=1
αj
N
X
k=1 |hˆu(tj),ΦkiV|2
subject to hΦi,ΦjiV=δij,1≤i, j ≤M.
93
CHAPTER 7. THE NONSTATIONARY CASE
Consider the linear mapping FN:RN→ VN, ek7→ ˆuk:= ˆu(tk), where
ekdenotes the k−th canonical basis vector of RNand ˆukare the snapshots.
Considering the following proposition, we can see that RNcorresponds with
Vand VNwith W. For a proof of the following singular value decomposition
(SVD), we refer for instance to [105].
Proposition 7.3. Let F:V→Wbe a linear operator, where Vand W
denote two finite-dimensional real Hilbert spaces with inner products h·,·iV
and h·,·iWand dimV =mand dimW =nwith m≥n. Then there exist
real numbers σ1≥σ2≥. . . ≥σn≥0and orthonormal bases {vk}n
k=1 of V
and {wk}n
k=1 of W, such that
F(vk) = σkwk, F∗(wk) = σkvk,
for k= 1,··· , n, where the adjoint operator F∗of Fis defined by the follow-
ing definition.
Definition 7.4. Let {V, h·,·iV}and {W, h·,·iW}be real Hilbert spaces and
F:V→Wa linear operator. We call F∗the adjoint operator of Fif
hw, FviW=hF∗w, viV
for all w∈Wand all v∈V.Fis called self-adjoint, if
F∗=F.
Then, we obtain with hv, wiRN=
N
P
j=1
αkvkwkfor all v∈RN
FN(v) =
N
X
j=1
αjhv, ejiRNFN(ej) =
N
X
j=1
αjhv, ejiRNˆuj.(7.2.2)
Assuming F∗
Nas the adjoint of FN, then follows
hFN(v),ΦiV=h
N
X
k=1
αkhv, ekiRNˆuk,ΦiV
=
N
X
k=1
αkhˆuk,ΦiVhek, viRN=
N
X
k=1
αkhˆuk,ΦiVvk
=h
hˆu1,ΦiV
.
.
.
hˆuN,ΦiV
, viRN
94
7.2. PROPER ORTHOGONAL DECOMPOSITION POD
for all Φ∈ VNand so, we can interpretate the adjoint operator as
F∗
NΦ =
hˆu1,ΦiV
.
.
.
hˆuN,ΦiV
(7.2.3)
for all Φ∈ VN. The idea is now to define RN:= FNF∗
Nand KN:=
F∗
NFN. Together with (7.2.2) and (7.2.3), we derive with hˆuk,FN(·)iV=
h
hˆuk,ˆu1iV
.
.
.
hˆuk,ˆuNiV
,·iRNand the following remark, see also [111], page 25,
RN=
N
X
j=1
αjhˆuj,·iVˆuj
KN=
hˆu1,FN(·)iV
.
.
.
hˆuN,FN(·)iV
.
Remark 7.5.
•The operator RNis bounded, self-adjoint and non-negative.
•By Hilbert-Schmidt theory exists an orthonormal basis {Φk}N
k=1 and
non-negative real numbers {λk}N
k=1 such that
RNΦi=λiΦiand λ1≥λ2≥. . . ≥λN.
By the Lagrangian theory it follows that the first-order optimality condi-
tion for the least-square problem (7.2.1) is
RNΦi=λiΦiand λ1≥λ2≥. . . ≥λN.(7.2.4)
For more details, see [106]. We are able to obtain the orthonormal basis
{Φk}N
k=1 by solving (7.2.4), the solution of (7.2.1).
Additionally, we have the following theorem to solve (7.2.1) by choosing
a fixed M.
Theorem 7.6. Let λ1≥λ2≥. . . ≥λN≥0be the non-negative eigenvalues
and {Φk}N
k=1 the associating eigenvectors of RN. Let MN, then {Φk}N
k=1
is orthonormal with rank Mand (7.2.1) satisfies
N
X
j=1
αjkˆu(tj)−
M
X
k=1hˆu(tj),ΦkiΦkk2
V=
N
X
j=M+1
λj.
95
CHAPTER 7. THE NONSTATIONARY CASE
For a proof, we refer to [111].
We see by Proposition 7.3 that the two orthonormal systems of the finite-
dimensional spaces can be transformed into the other one by a linear mapping
or its adjoint, if they are known.
Additionally, we find an orthonormal basis {vk}N
k=1 in RNsuch that for
k= 1, . . . , N
KN(vk) = λkvk.
Now, we are able to determine the optimal orthonormal basis {Φk}M
k=1 of
(7.2.1) in VNby the linear mapping FN
FN(vk) = pλkΦk,i.e. Φk=1
√λkFN(vk)
for a fixed M,k= 1, . . . , M.
In the following, we consider the problem of finding the so-called ’modes’
{ui(x)}M
i=0 so that the Galerkin approximations, which are defined with the
corresponding mean flows u0(x)=1/N PN
i=1 ˆui(x)by
u[M](x, t) :=
M
X
i=0
ai(t)ui(x),
with a0≡1and ai(t) := (u−u0, ui)Ω,minimizes the energy-related error
χu:= 1
N
N
X
i=1 kˆui(·)−u[M](·, ti)kL2(Ω)
compared to all other bases {wi(x)}M
i=1 and corresponding Galerkin approx-
imations, i.e.
χu≤χw.
For a homogeneous fluid and an incompressible flow, the flow velocities
u(x, t),having components uiin the xicoordinate direction, can be split-
ted into a mean part u(x)and a fluctuating part u0(x, t)using the so-called
Reynolds decomposition:
ui=ui+u0
i.
Now, u0(x)represents the the mean part u(x)and PM
i=1 ai(t)ui(x)without
u0(x)the fluctuation part u0(x, t)of the Reynolds decomposition.
This theory leads to the following algorithm to calculate the POD modes
and coefficients.
Algorithm 7.7.
96
7.2. PROPER ORTHOGONAL DECOMPOSITION POD
1. Compute the averaged (mean) flow
u0(x) := 1
N
N
X
i=1
ˆui(x).
We denote by ˜u=u−u0the fluctuation of ufrom this mean flow.
2. Compute the correlation matrix C∈RN×N
Ci,j =1
N(˜ui,˜uj).
3. Compute the eigenvalues λiand the set of normalized eigenvectors ~vi,
i= 1, . . . , N of Cwith :
C~vi=λi~vi.(7.2.5)
4. Compute the POD modes
ui:= 1
√Mλi
N
X
i=1
~vi˜ui.
5. Compute the Fourier coefficients
ai(tj) := (˜uj, ui)Ω.
For the theory in the infinite-dimensional case and the optimality of the
pod basis system, we refer to [106, 111].
Let us now present the POD modes for the Setting 7.1.
First, for the unactuated system, N= 567 snapshots ˆun
i(x) := ˆun(x, ti)
were determined at equidistant discrete times ti,i= 1,··· , N, covering
6convective time units. Analogously, Nsnapshots ˆua
i(x) := ˆua(x, ti),i=
1,··· , N, are computed for the actuated system by a URANS simulation with
aWilcox98 turbulence k-ω-model and a Reynolds number of 1.756 ·106.
We chose a fairly large actuation amplitude Bfor the actuated case to get
significant differences between the frequencies of the operating conditions.
The POD method for our Setting 7.1 yields the eigenvalues in Figure 7.2,
where the typical pairs of eigenvalues are demonstrated. In Figure 7.2, one
can also see that the first pair of modes contains the most energy. Due to
this reason, we consider only the first pair in the next chapter to introduce
the reduced-order model.
The mean flows is presented in Figure 7.3, the first mode in Figure 7.4
and the second mode is presented in Figure 7.5.
97
CHAPTER 7. THE NONSTATIONARY CASE
Figure 7.2: Eigenvalues of the natural (top) and the actuated (bottom) flow.
98
7.2. PROPER ORTHOGONAL DECOMPOSITION POD
Figure 7.3: Mean flow of the natural un
0= (un
0vn
0)T(un
0: top left, vn
0: bottom
left) and the actuated ua
0= (ua
0va
0)T(ua
0: top right, va
0: bottom right) case.
99
CHAPTER 7. THE NONSTATIONARY CASE
Figure 7.4: The first mode of the natural un
1= (un
1vn
1)T(un
1: top left, vn
1:
bottom left) and the actuated ua
1= (ua
1va
1)T(ua
1: top right, va
1: bottom right)
case.
100
7.2. PROPER ORTHOGONAL DECOMPOSITION POD
Figure 7.5: The second mode of the natural un
2= (un
2vn
2)T(un
2: top left,
vn
2: bottom left) and the actuated ua
2= (ua
2va
2)T(ua
2: top right, va
2: bottom
right) case.
101
CHAPTER 7. THE NONSTATIONARY CASE
102
Chapter 8
Reduced-order model (ROM)
As mentioned in the last chapter, we want to introduce a low-order model
describing the lift-increasing effect of high-frequency forcing.
In this chapter, we outline the approach of some engineers, see [66, 74], to
set up a reduced-order model without a complete and detailed mathematical
reflection. This procedure is based on many observations and is adapted to
the given problem related to Navier-Stokes equations and high-lift configu-
rations.
After a comprehensible explanation for the design of the reduced-order
model, adopted almost as it stands from [66, Section 3], we will present in
Section 8.1 a summary of the core statements of developing the ROM more
detailed, adopted almost as it stands from [66, Section 4and 5]. In section
8.2 and 8.3, we present our modifications on their dynamical system and a lift
formula based only on the Fourier coefficients. Numerical results to compare
this reduced-order model with the full turbulence model are considered at
the end of this chapter. The Sections 8.2-8.4 are published by John, Noack,
Schlegel, Tröltzsch and Wachsmuth in [57]. Based on this ROM, we will
establish in Chapter 9 a reduced optimization problem.
The dynamical system should reflect the following behavior of the un-
steady Reynolds-averaged Navier-Stokes (URANS) simulation:
(i) von K´arm´an vortex shedding without actuation: a vortex street is a phe-
nomenon of fluid mechanics for a repeating pattern of swirling vortices
behind a bluff body caused by the unsteady separation of fluid flow. It
is named after the engineer, Theodore von K´arm´an (1881-1963),
(ii) lock-in shear-layer shedding under high-frequency forcing: a shear-layer
is the transition region between two parallel fluid flows and a shear-
layer shedding means that a boundary-layer fixed to a body separates
from the body surface,
103
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
(iii) a transient behavior from the natural case (i)to the actuated one under
forcing (ii): this means that the system should describe the transition
from the natural flow to the actuated flow and
(iv) a transient behavior from the actuated case (ii)to the natural case (i)
when forcing is turned off.
Due to the periodic actuation, we have to consider oscillatory flows, which
are characterized by an amplitude Aand a phase α, i.e. the argument of
sinusoidal functions, if they are linear in time. We can consider them as
polar coordinates of the phase-space (a1, a2) = A[cos α, sin α].
First, we are searching for a system describing the natural flow. Noack [74]
described the self-amplified amplitude-limited behavior of vortex shedding
by the following Landau-equation (8.0.1). Let therefore, the superscript n
stand for the unactuated natural case, σnfor the positive growth rate, σn,n
for the positive Landau constant and An=pa2
1+a2
2for the amplitude. For
simplicity, we assume the frequency ωnas a constant.
˙a1= ˜σna1−ωna2
˙a2=ωna1+ ˜σna2
˜σn=σn−σn,n(An)2,
(8.0.1)
The superscript astands in the following for the actuated case.
The shear-layer dynamics is stimulated by high-frequency forcing g(t) =
Bcos(β)with an amplitude B, a phase βand the frequency ˙
β=ωa.The
shear layer denotes in fluid mechanics the transition area between two parallel
streams with different velocities in contrast to wall-bounded boundary layers.
The phase difference of the actuation with respect to the oscillation of the
flow is denoted by θ. That means that the oscillation flow has the phase
θ+β. The behavior is most easily modeled by a linear damped oscillator
with a periodic forcing at the eigenfrequency. Here is σaa negative growth
rate and g3, g4∈Rare parameters to describe the gain of actuation. In
contrast to the natural case(8.0.1), we use in the actuated case the indices 3
and 4instead of 1and 2.Thus, the actuated flow is described by the system:
˙a3=σaa3−ωaa4+g3Bcos(θ+β),
˙a4=ωaa3+σaa4+g4Bsin(θ+β).(8.0.2)
In reality, every flow ˜uwith actuation ˜gconsists of a superposition of
several frequencies ˜ω1,...,˜ωN,N∈N,such that its energy E(˜g)is the sum
of energies associated to the frequencies Ei(˜g),i= 1, . . . , N :
E(˜g) =
N
X
i=1
Ei(˜g).
104
Notice that we assume in both the natural flow, gn= 0,and the actuated flow,
ga=Bacos(ωat),a constant frequency ωnand ωa, respectively, such that the
associated energies E(gn)and E(ga)are related to the only frequencies ωn
and ωa.Thus, we characterize the energies in the natural and the actuated
flow by the squares of the amplitudes (An)2and (Aa)2,respectively:
E(gn) := (An)2and E(ga) := (Aa)2.
Let us analogously denote by (Ai
˜g)2the energy of the flow with actuation ˜g
associated to the frequency ωi:
(Ai
˜g)2:= Ei(˜g).(8.0.3)
Now, we comprise both oscillations in a four-dimensional phase space
[a1a2a3a4]. With B= 0, i.e. a3=a4= 0, this system describes the natural
flow, according to (i). To obey (ii), the oscillation at the natural frequency
has to vanish with increasing actuation amplitude B > 0. To achieve this,
we decrease the growth rate of the natural case σnwith the help of the
growth of the high-frequency amplitude Aa=pa2
3+a2
4. Analogously to the
damping term −σn,n(An)2in the Landau system, we add additionally the
term −σn,a(Aa)2and we get
˜σn=σn−σn,n(An)2−σn,a(Aa)2.
We see that the energy in the natural case decreases with increasing
energy in the actuated case and vice versa. This guarantees a1=a2= 0 at
the forced state, according to (ii). We substitude the terms g3Bcos(θ+β)
and g4Bsin(θ+β)in (8.0.2) by g31g+g32 ˙gand g41g+g42 ˙g, respectively, to
guarantee more flexibility by calibrating this system to original results, see
Section 8.2.
Thus, we introduced a low-order dynamical system of two coupled os-
cillators describing the observed behavior of the natural and the actuated
flows as well as a transient behavior between them, according to the desired
properties (i)−(iv):
˙a1= ˜σna1−ωna2
˙a2=ωna1+ ˜σna2
˙a3= ˜σaa3−ωaa4+g31g+g32 ˙g
˙a4=ωaa3+ ˜σaa4+g41g+g42 ˙g
˜σn=σn−β1(An)2−β2(Aa)2
˜σa=σa
(8.0.4)
105
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
with An=pa2
1+a2
2,Aa=pa2
3+a2
4. The ai’s, i= 1,...,4can be inter-
preted as the coefficients to the in Chapter 7 calculated POD modes. With
the steady base flow u0,this leads to
u(x, t)≈u0(x) +
4
X
i=1
ai(t)ui(x).
If uis a velocity field without actuation, a3and a4are almost equal to
zero. Thus, uis approximated by the modes of the natural case u(x, t)≈
u0(x) + a1(t)u1(x) + a2(t)u2(x).Otherwise, uis described by modes of the
actuated case and a1and a2are near to zero. One can see in the system
(8.0.4) that more actuation, e.g. for instance a higher actuation amplitude
B, leads to a higher weighting of the modes of the actuated case in contrast
to the ones of the natural case.
Here, we see a low-order model with the control function g(t). This model
replace the Navier-Stokes equations in our reduced optimization problem, see
Chapter 9.
After a more detailed discussion for the reduced-order model (ROM) in
Section 8.1, adopted from [66], we explain in the following sections our mod-
ifications on this model (Section 8.2), introduce a lift-formula only based on
the mode coefficients {ai}4
i=1 (Section 8.3) and give a numerical example to
demonstrate that the ROM reproduces the nonlinear behavior of the system
sufficiently well for our optimization ansatz (Section 8.4).
8.1 A generalized model
This section, where we consider the structure of the dynamical system (8.0.4)
more detailed, is based on [66, Section 4 and 5].
8.1.1 Mean-field theory
We consider a computational domain Ω⊂R2with x= (x, y)∈Ω.The
x-axis is aligned with the flow and the y-axis with the orthogonal direction.
The velocity field is denoted by u= (u, v),where the components uand v
are aligned with the x- and y-direction. This model will demonstrate the
role of mean-field dynamics in stabilizing an attractor and as a commitment
between the actuated and the natural (unactuated) case. An attractor is a
set towards which a dynamical system evolves over time1. For the mean-field
1Wikipedia
106
8.1. A GENERALIZED MODEL
theory, we refer to [74, 75, 94, 95]. In addition to the L2-scalar product of
vector fields
(f, g)Ω=ZΩ
f·gdx,
let us define inner product for matrix-valued fields by
(A, B)Ω:= Z
Ω
A:Bdx, with A:B:=
2
X
i,j=1
AijBji
and the instantaneous kinetic energy of a velocity field uis given by
K:= kuk2
Ω
2.
We consider the incompressible
div u= 0
non-stationary Navier-Stokes equation
ut−ν∆u+ (u·∇)u+∇p= 0 (8.1.1)
and the unsteady boundary condition
u=g
with boundary actuation g. Additionally to the boundary actuation, we
consider a time periodic and space dependent volume force ga. We denote
the so-called ensemble average by uwith its approximation
u(t) := 1
TZT/2
−T/2
u(t+τ)dτ,
where T > 0is a set length of a time window. Next, we are formulating some
assumptions. The first one is based on observed phenomenology.
(A.1) (A generalized Krylov-Bogoliubov ansatz) The velocity field u
is dominated by the sum of a slowly varying base flow uband two
oscillatory components which are nearly pure harmonics at the natural
unand the actuation uafrequency. Other temporal harmonics are
considered as negligible. Thus, we obtain
u(x, t) = ub(x, t) + un(x, t) + ua(x, t),(8.1.2)
where ubsatisfies the steady, inhomogeneous boundary condition g(t),
unthe homogenized version and uaaccounts for the residual to the
unsteady boundary condition g(t)−g(t).
107
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
Due to the fact that the base flow is almost time-independent, we assume
ub=ub.
Furthermore, we recognize that the ensemble averages of the fluctuations un
and uavanishes
un= 0, ua= 0
where ub,unand uaare averaged over the associated time intervals.
The ansatz of Dušek et. al. [38] was to establish a small parameter 1
and slowly varying amplitude functions ub
0,un
iand ua
i,i= 1,2, such that
ub(x, t) = ub
0(x, t),
un(x, t) = un
1(x, t) cos(Ωnt) + un
2(x, t) sin(Ωnt),
ua(x, t) = ua
1(x, t) cos(Ωat) + ua
2(x, t) sin(Ωat),
expresses the assumed slow variation of the mean flow, the oscillation ampli-
tudes, the frequencies and the phase shifts.
This ansatz implies that time derivatives of the amplitude functions are
of order O(), which we want to neglect.
The second assumption is called by engineers ’a non-commensurability
ansatz’:
(A.2) (A non-commensurability ansatz) There is no direct interaction
between unand uathrough the nonlinear term (u·∇)u, i.e. (un·∇)ua=
(ua·∇)un= 0.
This assumption is based on the numerically observed fact that the activity
regions of these fluctuations rarely overlap. So, on each of the two attractors,
we neglect fluctuations in the other frequency.
Substituting the Assumption (A.1) into the Navier-Stokes equations (8.1.1)
and re-arranging the terms by the zeroth and the first harmonics at frequen-
cies Stnand Sta, respectively, leads to
0 = −ν∆ub+ (ub·∇)ub+ (un·∇)un+ (ua·∇)ua+∇pb,(8.1.3)
un
t=ν∆un−(un·∇)ub−(ub·∇)un−∇pn,(8.1.4)
ua
t=ν∆ua−(ua·∇)ub−(ub·∇)ua−∇pa+ga.(8.1.5)
The temporal behavior of the terms (uk· ∇)uk,k=b, n, a, is specified by
the 0th and second harmonics of the frequencies Stnand Staand they are
eliminated in (8.1.4) and (8.1.5) by Assumption (A.1). The mean-field model
in the next subsection is based on the system (8.1.3)-(8.1.5).
108
8.1. A GENERALIZED MODEL
For a homogeneous fluid and an incompressible flow, the flow velocities
u(x, t),having components uiin the xicoordinate direction, can be splitted
into a mean part u(x)and a fluctuating part u0(x, t)using the so-called
Reynolds decomposition:
ui=ui+u0
i.
The Reynolds stress tensor τ0is then defined by τ0
ij =u0
iu0
jand describes the
degree of nonlinearity. The term ∇τ0denotes the force pushing the mean
flow away from the steady flow. If the Reynolds stress would be zero, then
the mean flow would coincide with the steady flow.
The next and last assumption guarantees a linear relation between the
Reynolds stresses τ0and the mean-field correction term uh. Let therefore, us
be the associating solution of the steady Navier-Stokes equations
−∆us+ (us·∇)us+∇p= 0 in Ω
div us= 0 in Ω
us=g(t)on Γ.
(8.1.6)
(A.3) (Linearized Reynolds equation) Let us assume that (8.1.3) is lin-
earizable at the steady solution usand
ub=us+uh.(8.1.7)
The linearized Reynolds equation for the mean-field correction uhis
obtained by substituting (8.1.7) in (8.1.3)
0 = −ν∆(us+uh) + ((us+uh)·∇)(us+uh)
+ (un·∇)un+ (ua·∇)ua+∇(ps+ph)
subtracting the steady Navier-Stokes equation
0 = −ν∆uh+ (us·∇)uh+ (uh·∇)us+∇ph
−ν∆us+ (us·∇)us+∇ps−(−ν∆us+ (us·∇)us+∇ps)
+ (un·∇)un+ (ua·∇)ua+ (uh·∇)uh
and neglecting quadratic terms in uh:
0 = −ν∆uh+ (us·∇)uh+ (uh·∇)us
+(un·∇)un+ (ua·∇)ua+∇ph.(8.1.8)
109
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
8.1.2 Mean-field Galerkin model
Let us now transform the mean-field model over to the least-order Galerkin
model. For the theory of the Galerkin method, we refer to to Fletcher [41],
Holmes et. al. [56] and Ladyzhenskaya [65]. Let u0denote the steady base
flow, then the Galerkin method is based on the Galerkin expansion
u(x, t) =
M
X
i=0
ai(t)ui(x).(8.1.9)
The time dependency is described by the Fourier coefficients ai, i = 1, . . . , M
with a0= 1.We describe the Galerkin approximation in the section Least-
order Galerkin approximation after next.
Galerkin method
Replacing the flow uin the Navier-Stokes equations (8.1.1) with ga= 0 by the
Galerkin expansion (8.1.9) and projecting them onto the subspace spanned
by the expansion modes:
(ut−ν∆u+ (u·∇)u+∇p, ui)Ω= 0 for i= 1, . . . , M,
we obtain with Γ := ∂Ω
1.
∂
∂t M
X
j=0
aj(t)uj(x)!, ui(x)!Ω
= M
X
j=0
∂aj
∂t uj, ui!Ω
=
M
X
j=1 ∂aj
∂t uj, uiΩ
=
M
X
j=1
∂aj
∂t (uj, ui)Ω=
M
X
j=1
∂aj
∂t δij =∂ai
∂t ,
2.
ν ∆ M
X
j=0
aj(t)uj(x)!, ui(x)!Ω
=ν
M
X
j=0
(∆(ajuj), ui)Ω
=ν
M
X
j=0
aj(∆uj, ui)Ω
| {z }
=:lij
=ν
M
X
j=0
lijaj,
110
8.1. A GENERALIZED MODEL
3. M
X
j=0
aj(t)uj(x)!·∇! M
X
k=0
ak(t)uk(x)!, ui(x)!Ω
=
M
X
j=0
M
X
k=0
(((ajuj)·∇)(akuk), ui)Ω
=
M
X
j=0
M
X
k=0
ajak((uj·∇)uk, ui)Ω
| {z }
=:qu
ijk
=
M
X
j,k=0
ajakqu
ijk,
4. By Noack et. al. [76], we obtain a pressure expansion ∇p(x, t) =
p[0...M](x, t) =
M
P
j=0
M
P
k=0
pjk(x)aj(t)ak(t). This leads to
(∇p(x, t), ui(x))Ω=p[0...M](x, t), ui(x)Γ
= M
X
j=0
M
X
k=0
pjk(x)aj(t)ak(t), ui(x)!Γ
=
M
X
j=0
M
X
k=0
ajak(pjk, ui)Γ
| {z }
=:qp
ijk
=
M
X
j,k=0
qp
ijkajak.
Summarized, this leads with qijk := qu
ijk +qp
ijk to a simplified ordinary differ-
ential equation system
∂ai
∂t =ν
M
X
j=0
lijaj+
M
X
j,k=0
qijkajak
i= 1, ..., M, to define the Fourier coefficients ai(t). Together with
(M(a, a))i:=
M
X
j,k=1
qijkajak,(F(a))i:=
M
X
j=1
(νlij +qi0j+qij0)aj,
(C)i:= νli0+qi00
and a0= 1,we obtain the system
∂a
∂t =M(a, a) + F(a) + C.(8.1.10)
We have the following difficulty: by a standard POD method, boundary
values is prescribed by the dynamical system and can not be installed as a
111
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
actuation command. The problem is that boundary actuation is generally
not derivable from the Galerkin projection, because the Galerkin projection
is composed to ignore boundary perturbations defined over a set of measure
zero. But it is possible to simulate boundary effects by additional actuation
modes, see for instance [16, 49, 77].
Instead, we decided to add an control term to the system (8.1.10) that is
state dependent and includes the control command and its time derivative.
We consider an oscillatory control with slowly varying periodic characteris-
tics.
Then the influence on the flow can be modelled by an actuation term
B~g(t)with actuation command ~g(t)and a matrix B, see [85, 86, 92]. As
mentioned in (7.0.3), we consider a periodic actuation g(t) = Bcos(β(t)),
fulfilling ∂β(t)/∂t =−Ωawith an actuation amplitude Band a phase shift
of β(t)−Ωat. The acceleration is ∂g(t)/∂t =−ΩaBsin(β(t)).
We combine the actuation command and its derivative with respect to t
to
~g = [g, ˙g]T=B[cos(β),−Ωasin(β)]T
and obtain the Galerkin system with actuation
∂a
∂t =M(a, a) + F(a) + C+B~g. (8.1.11)
The term B~g(t)replaces the boundary actuator in the Navier-Stokes equa-
tions, because we have no boundary terms in the dynamical system (8.1.11).
Considering the system (8.0.4), we have to define Bby
B=
0 0
0 0
0 0
0 0
g31 g32
g41 g42
.
Let us consider in the next section the Galerkin approximation for our
problem.
Least-order Galerkin approximation
The least-order Galerkin approximation is based on Assumption (A.1); we
are interested in modes resolving un,uaand ubin (8.1.2). Let therefore un
i
and ua
i,i= 1,2, be the two dominant POD modes of the natural and the
actuated attractors. The modes un
1and un
2are orthonormal and un
1and un
2
are orthonormal by construction.
112
8.1. A GENERALIZED MODEL
Considering the computed POD mode pairs, see Figure 7.3-7.5, we rec-
ognize that the actuated modes ua
1and ua
2have their main fluctuations over
and near the airfoil, whereas the fluctuations of the natural POD modes un
1
and un
2are concentrated further downstream. This fact and the different
wavelengths shows that POD mode pairs are nearly orthogonal. This ob-
servations supported by very small values (un
1, ua
1)Ω,(un
1, ua
2)Ω,(un
2, ua
1)Ωand
(un
2, ua
2)Ω.
Thus, we merge them into an orthonormal basis, using the Gram-Schmidt
normalization and obtain (u1, u2, u3, u4)associated to (un
1, un
2, ua
1, ua
2)and we
shall maintain the approximation of the fluctuations unand uaby
un=
2
X
i=1
ai(t)ui(x),and ua=
4
X
i=3
ai(t)ui(x).
Now, we are looking for a representation for ub. Let usbe the steady
solution of (8.1.6), un
0is the mean flow of the natural attractor and x∝y
means that xis proportional to y. Then, following [74], the effect of the
Reynolds stress due to the natural oscillations is described by a shift-mode
un
∆∝un
0−us,. Analogously, we define ua
∆.
Assume that u5and u6are obtained by a Gram-Schmidt orthogonal-
ization from un
∆and ua
∆, removing any projections over u1,··· , u4. Thus,
u1, . . . , u6are orthonormal. In [74], they approximate the time-varying base
flow ub(x, t) = us(x) + uh(x, t), see Assumption (A.3), with the two shift-
modes u5and u6corresponding to the two attractors of the natural and the
actuated case
ub(x, t) = us(x) + uh(x, t) = us(x) + a5(t)u5(x) + a6(t)u6(x),(8.1.12)
This means that the fluctuations u1, . . . , u4are negligible to approximate the
mean-field correction uhand the base flow ub.
A transient describes the crossing behavior into another attractor without
reaching him. Then, the mean flows un
0and ua
0are approximated by the
associated initial and final values of the base flow ubtrajectory in transients
connecting the two attractors and we obtain with u∆:= (un
0−ua
0)/kun
0−ua
0kΩ
the approximation
ub(x, t) = un
0(x) + a∆(t)u∆(x).(8.1.13)
See Figure 8.1 for a visual description of the relation between the steady
solution usand the mean flows un
0and ua
0.
113
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
Figure 8.1: The relation between the steady solution usand the mean flows
un
0and ua
0, see [66] Figure 5.
Considering (8.1.13) instead of (8.1.12), the velocity field of the URANS
data was approximated by
u(x, t) = un
0(x) + a∆(t)u∆(x) +
4
X
i=1
ai(t)ui(x).(8.1.14)
Thus, we do not have to extract a steady solution us.
In (8.1.14), we approximated the system of the natural flow, the actuated
flow and the states between them by the two mean flows un
0and ua
0, the
associated shift mode u∆with Fourier coefficient a∆, as transition between
them, the two modes u1, u2of the natural and the two modes u3, u4of the
actuated flow, with associated Fourier coefficients a1, ..., a4.
However, in the next section, we use (8.1.12) to obtain the equivalent
Galerkin expansion
u(x, t) = us(x) +
6
X
i=1
ai(t)ui(x).(8.1.15)
This approximation will be used in the Galerkin system (8.1.16)-(8.1.18).
The main advantages of (8.1.15) in contrast to (8.1.14) appear by cali-
brating the Galerkin system parameter with empirical data, see [66] for more
114
8.1. A GENERALIZED MODEL
details. The advantage of (8.1.14) was that we do not have to calculate the
steady solution us.
Least-order Galerkin system
Inserting the equation (8.1.15) into the mean-field Navier-Stokes equations
(8.1.4), (8.1.5) and (8.1.8), followed by the Galerkin projection of these equa-
tions on the expansion modes lead to the least-order Galerkin system, con-
sisting of the Fourier coefficients ai. Subsequently, we have to enforce a
Galerkin projection of these equations on the expansion modes {ui}6
i=1.
The Galerkin system, associated to the projection of (8.1.4), is
∂ai
∂t
|{z}
b=∂tun
=
2
X
j=1
6
X
k=5
qijkajak
| {z }
b=−(un·∇)ub
+
6
X
j=5
2
X
k=1
qijkajak
| {z }
b=−(ub·∇)un
+ν
2
X
j=1
lijaj
| {z }
b=ν∆un
+ 0
|{z}
b=∇pn
,
i∈ {1,2}.For a proof, see the derivation of(8.1.10). The other equations
follow analogously.
Let eibe the unit-vector of RM,then we get with
an=a1e1+a2+e2, aa=a3e3+a4e4,and ab=a5e5+a6e6
the full Galerkin system by
∂an
∂t =M(ab, an) + M(an, ab) + F(an),(8.1.16)
∂aa
∂t =M(ab, aa) + M(aa, ab) + F(aa) + B~g, (8.1.17)
0 = M(an, an) + M(aa, aa) + F(ab).(8.1.18)
In [66, Appendix B], it is shown that the Galerkin system has the following
equivalent form
∂
∂t
a1
a2
a3
a4
=
˜σn−˜ωn0 0
˜ωn˜σn0 0
0 0 ˜σa−˜ωa
0 0 ˜ωa˜σa
a1
a2
a3
a4
+
0 0
0 0
κ−η
η κ
~g, (8.1.19)
with ˜σn=σn−σn,n(An)2−σn,a(Aa)2,
˜ωn=ωn+ωn,n(An)2+ωn,a(Aa)2,
˜σa=σa−σa,n(An)2−σa,a(Aa)2,
˜ωa=ωa+ωa,n(An)2+ωa,a(Aa)2,
(8.1.20)
115
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
where An=kank,Aa=kaakare the respective oscillation amplitudes. The
Assumption (A.3) avoids a Taylor series in (An)2and (Aa)2.
The parameter σn,n describes the decreasing growth rate of the natural
attractor with increasing energy in the natural attractor and σn,a the de-
creasing growth rate of the natural attractor with increasing energy in the
actuated attractor. The parameter σa,n and σa,a can be interpreted analo-
gously.
The parameters ωn,n, ωn,a, ωa,n and ωa,a describe the changes of the amplitude-
dependent frequency. The equation (8.1.18) leads to a linear dependence of ab
on anand aa.Consequences of this dependency are that the system (8.1.19)
only consists of a1, . . . , a4and the equations (8.1.20), for more details see [66,
Appendix B].
Let us have a closer look at (8.1.19) and (8.1.20) to identify some pa-
rameters. A first ansatz is that we assume for reasons of simplifications
that the oscillation frequencies are independent of Anand Aaand a constant
frequency for the natural and the actuated flow, which yields
ωn,n =ωn,a =ωa,a =ωa,n = 0
and consequently
˜ωn=ωn= Ωn,
˜ωa=ωa= Ωa,
see [66, Section 5.4] for more details. In [66, Section 5.5] is the simplifi-
cation
˜σa=σa,
described that means σa,n =σa,a = 0.
Finally, this leads to the following dynamical system
˙a1= ˜σna1−Ωna2
˙a2= Ωna1+ ˜σna2
˙a3= ˜σaa3−Ωaa4+κg −η˙g
˙a4= Ωaa3+ ˜σaa4+ηg +κ˙g
˜σn=σn−σn,n(An)2−σn,a(Aa)2
˜σa=σa
(8.1.21)
and one has to calibrate the remaining parameters to the numerical data,
i.e. the preliminarily calculated snapshots.
Similar to the least-order Galerkin model (8.1.19) and (8.1.20), we want
to establish a new reduced-order model (ROM) calibrated with our numerical
data. In the next section, we present our modifications.
116
8.2. MODIFICATIONS ON THE REDUCED-ORDER MODEL
8.2 Modifications on the reduced-order model
In this section, we want to demonstrate our modifications and calibrations
on the reduced-order model to handle the optimization problem (PN)in the
next chapter.
Similarly to POD, see Section 7.2, all snapshots are processed. For this
purpose, we consider only the velocity field u= (u, v)in a certain reference
domain, where the actuation has the main influence on velocity and lift, see
Figure 7.1. Therefore, the snapshot velocity data are weighted by the size of
their area. We select the first two POD modes of the actuated (ua
1, ua
2)and
non-actuated (un
1, un
2)system, carrying the highest energy.
But in contrast to the ansatz in the last Section 8.1; we filter the POD
modes and the associated mode coefficients by eliminating certain fluctua-
tions to emphasize the dominant harmonic structure, see the following sub-
section for more details.
8.2.1 Filtering of the POD modes and coefficients
The snapshots ˆun
iand ˆua
iare given on a time interval [0, T]. Due to the fact
that the natural flow and the actuated flow have different wavelengths, we
search for the maximal ka, kn∈Nsuch that the times
Tn= 2πkn/Ωnand Ta= 2πka/Ωa
fulfill Tn< T and Ta< T. Therefore, the time Thas to be big enough such
that T > 2π/Ωnand T > 2π/Ωa.
Let (an
1(t), an
2(t)) and (aa
1(t), aa
2(t)) are the first POD mode coefficient
pairs of the natural and the actuated case, respectively. We recognize that
they are very similar to trigonometric functions or oscillations and so, we
want to approximate them as well as possible by oscillations. Therefore, we
calculate the phases ϕn(t),ϕa(t)and radii ˜rn(t),˜ra(t)by
an
1(t) + ian
2(t) = ˜rn(t)eiϕn(t),
aa
1(t) + iaa
2(t) = ˜ra(t)eiϕa(t).
To extract the dominant harmonic oscillations from these POD coefficients,
we smooth in (8.2.1) perturbations of both, the radii ˜rn(t),˜ra(t)and the
phases ϕn(t), ϕa(t), in anticipation of a small deformation of the correspond-
ing modes. This holds, because the POD method does not extract pure
frequencies and radii, but ’deformed’ modes with the highest energy.
With the average values
rn= ˜rn(t) := 1/TnRTn/2
−Tn/2˜rn(t)dt, ra= ˜ra(t),
ωn=∂tϕn(t), ωa=∂tϕa(t),(8.2.1)
117
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
we approximate our filtered coefficients
a1(t) + ia2(t) = rneiωnt,
a3(t) + ia4(t) = raeiωat.
Let un
0and ua
0are the mean flows
un
0(x) = 1/N
N
X
i=1
ˆun
i(x)
in the unactuated and
ua
0(x) = 1/N
N
X
i=1
ˆua
i(x)
in the actuated case.
Then, the associated filtered modes are determined by the Fourier ansatz
ui(x) = ˆun(x, ·)−un
0,1
rnai(·)Tn
:= 1
Tn
Tn/2
Z
−Tn/2
(ˆun(x, t)−un
0(x)) 1
rnai(t)dt,
for i= 1,2and
ui(x) = ˆua(x, ·)−ua
0,1
raai(·)Ta
=1
Ta
Ta/2
Z
−Ta/2
(ˆua(x, t)−ua
0(x)) 1
raai(t)dt,
for i= 3,4.The coefficients an
1and an
2are orthonormal and aa
1and aa
2are
orthonormal by construction. But the coefficients an
iare not necessarily or-
thonormal to aa
ifor i= 1,2.Thus, we finally orthonormalize these frequence-
filtered modes ui,i= 1,...,4by Gram-Schmidt and denote the orthonormal
modes for simplicity by ui,i= 1,...,4.
The Figures 8.6 and 8.7 show the associated filtered mode coefficients a1,
i= 1,...,4in contrast to the original ones.
In the Figures 8.2 and 8.3, we present the filtered modes u1and u2of the
natural flow and in 8.4 and 8.5, we present the filtered modes u3and u4of
the actuated flow in comparison to the original POD modes. We see that the
filtered modes of the actuated case u3and u4have their main fluctuations
over and near the airfoil, whereas the fluctuations of the filtered modes of
the natural case u1and u2are concentrated further downstream. Because we
approximated the mode coefficients a1,i= 1,...,4to almost pure trigono-
metric functions, we see that some original modes are smoother than the
filtered ones.
118
8.2. MODIFICATIONS ON THE REDUCED-ORDER MODEL
Figure 8.2: The first filtered mode u1= (u1v1)T(u1: top left, v1: bottom
left) and the original first mode un
1= (un
1vn
1)(un
1: top right, vn
1: bottom
right) of the natural case.
119
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
Figure 8.3: The second filtered mode u2= (u2v2)T(u2: top left, v2: bottom
left) and the original first mode un
2= (un
2vn
2)(un
2: top right, vn
2: bottom
right) of the natural case.
120
8.2. MODIFICATIONS ON THE REDUCED-ORDER MODEL
Figure 8.4: The first filtered mode u3= (u3v3)T(u3: top left, v3: bottom
left) and the original first mode ua
1= (ua
1va
1)(ua
1: top right, va
1: bottom
right) of the actuated case.
121
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
Figure 8.5: The second filtered mode u4= (u4v4)T(u4: top left, v4: bottom
left) and the original first mode ua
2= (ua
2va
2)(ua
2: top right, va
2: bottom
right) of the actuated case.
122
8.2. MODIFICATIONS ON THE REDUCED-ORDER MODEL
Figure 8.6: The filtered mode coefficients (natural case a1: top left, actuated
case a3: bottom left) and the original mode coefficients (natural case an
1: top
right, actuated case aa
1: bottom right) of the associated first mode over the
snapshots.
123
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
Figure 8.7: The filtered mode coefficients (natural case a2: top left, actuated
case a4: bottom left) and the original mode coefficients (natural case an
2: top
right, actuated case aa
2: bottom right) of the associated second mode over
the snapshots.
8.2.2 Parameter calibration
The filtered modes (u1,...u4)contain significant information gained from
the URANS solution by the k-ω-model. Considering the dynamical system
(8.1.21), we redefined β1:= σn,n and β2:= σn,a and decided to replace κand
νby g31, g42 and g32, g41, respectively, to have more degrees of freedom to
calibrate the reduced-order model to the original URANS data.
124
8.2. MODIFICATIONS ON THE REDUCED-ORDER MODEL
This leads to the following reduced-order model:
˙a1= ˜σna1−ωna2
˙a2=ωna1+ ˜σna2
˙a3= ˜σaa3−ωaa4+g31g+g32 ˙g
˙a4=ωaa3+ ˜σaa4+g41g+g42 ˙g
˜σn=σn−β1(An)2−β2(Aa)2
˜σa=σa,
(8.2.2)
with An=pa2
1+a2
2,Aa=pa2
3+a2
4and g=Bcos(ωat), where Bis the
amplitude of the actuation signal and ωais the associated angular frequency.
Our snapshots contain no data or information about the transient behav-
ior, because the snapshots are taken while the natural flow and the actuated
one, respectively. Because of this reason, we have to select the amplification
rates σn,σaas follows by empirical values:
•σn= 0.15 is an empirical value, if the cord-length of the wing is 1.
Because the flap is the active part of the configuration, we choose σn=
0.15U∞
cfl .
•σa=−1
Tcon where Tcon is the time that one vortex needs to pass the
flap-length cfl. We read off this value from the snapshots.
Now, we want to calibrate the parameters β1and β2of the system (8.2.2).
They are determining the growth rates of both oscillations. That means that
they are responsible for the increasing or decreasing rate of the energies in
both oscillations with decreasing or increasing energy in the other one. If
the fluid flow is in the unactuated state, then no energy should be in the
coefficients a3, a4, i.e. a3=a4= 0, hence Aa= 0. Moreover, we require that
˜σn=σn−β1(An)2= 0
holds for unactuated flow dynamics. This expresses the fact that there is no
additional energy contribution to the natural oscillatory behavior of a1, a2.
Thus β1can be determined by
β1=σn1
(An)2=σn1
(rn)2.
There are many possibilities to calibrate β2. For instance, we can determine
β2analogously to β1by
β2=σn1
(Aa)2=σn1
(ra)2.
125
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
The problem of this ansatz is that the energy belonging to the natural fre-
quency An
ga, see (8.0.3), does not vanish completely in the actuated case.
That means An
ga>0.
Let us assume that the energy in the whole system is constant over all
actuation amplitudes B:
σn=β1(An)2+β2(Aa)2.(8.2.3)
Then, another possibility is to determine β2by
β2=σn−β1(An)2
(Aa)2.
We decided to follow an alternative approach. Let us consider for example
the actuation amplitude ˜
Bwith associated flow ˜u. Till now, we assumed a
constant frequency in both the natural and the actuated flow. To determine
the parameters β1and β2,we ignore this. Instead, we assume that every flow
consists of a combination of the natural and the actuation frequency. Thus,
the energy in this flow ˜uis the sum of the energies associated to the two
frequencies. Then, we denote analogously to (8.0.3) by (An
˜
B)2and (Aa
˜
B)2the
energies associated to the natural and the actuation frequency.
Let Babe the actuation amplitude related to the actuated flow ua. Con-
sidering the to the natural and the actuated flow associated energies (An
0)2,
(Aa
0)2,(An
Ba)2and (Aa
Ba)2,we are able to draw the energies over the actuation
amplitude B, see Figure 8.8.
Considering Figure 8.8, we only obtain a linear relation without con-
straints between the energies of the natural ωnand the actuation frequency
ωa,respectively, over the actuation amplitude B. But, this does not corre-
spond to reality. Therefore, we need an additional set of snapshots with a
small actuation amplitude B∗to identify the parameters β1and β2in the
system (8.2.2). We select B∗such that the associated energy belonging to
the natural frequency (An
B∗)2should be significant greater than zero. For
this actuation amplitude, we compute the filtered coefficients {a∗
i}4
i=1 and
the associated energies (An
B∗)2and (Aa
B∗)2.
We draw the results in Figure 8.9 and assume that (An
Ba)2is the lower
bound for the energy of the natural frequency ωnand (Aa)2is the upper
bound for the energy associated with the actuation frequency ωa.
Then, we determine the parameters β1and β2by assuming that the energy
in the whole system remains constant over all amplitudes B(8.2.3) by the
system
β1(An
0)2+β2(Aa
0)2=σn,
β1(An
B∗)2+β2(Aa
B∗)2=σn.(8.2.4)
126
8.2. MODIFICATIONS ON THE REDUCED-ORDER MODEL
Figure 8.8: Energies (An)2(continuous lines) and (Aa)2(dashed lines) over
the actuation amplitude B. Here, we calculated the energies only for B= 0
and the actuation amplitude Ba.
The reduced-order model (8.2.2) contains additionally free parameters
g31,g32,g41 and g42 to calibrate the selected actuation to the dynamical
system. Recall therefore that the actuation gand its derivative ˙gare g=
Bcos(ωat)and ˙g=−Bωasin(ωat),where the actuation amplitude Bis our
optimization variable.
Therefore, we multiply the third and fourth equation of (8.2.2) by gand
˙g, respectively, and integrate over [0, Ta]. This eliminates g32,g42 and g31,
g41, respectively. For instance
( ˙a3, g)Ta=σa(a3, g)Ta−ωa(a4, g)Ta+g31(g, g)Ta
leads to
g31 =( ˙a3, g)Ta−σa(a3, g)Ta+ωa(a4, g)Ta/(g, g)Ta.(8.2.5)
Note that (˙g, g)Tavanishes in the long term average.
An example of the phase portraits for the coefficients a1,··· , a4of the
system (8.2.2) is presented in Figure 8.10. This figure presents the solution of
(8.2.2) starting on the natural attractor with actuation. We see in this figure
on the left side in the phase portrait of the first oscillator (a1, a2), describing
127
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
Figure 8.9: Energies (An)2(continuous lines) and (Aa)2(dashed lines) over
the actuation amplitude B. Here, we calculated the energies for B= 0, B∗
and the actuation amplitude Ba.
the natural flow, that the energy vanishes and transfers over to the phase
portrait of the second oscillator (a3, a4),describing the actuated flow.
Figure 8.10: Phase portraits of (a1, a2)(left) and (a3, a4)(right) of the system
(8.2.2) with full actuation, starting with (a1, a2)on the natural limit cycle
and (a3, a4) = (0,0) with actuation.
128
8.3. COMPUTATION OF LIFT
The next Figure 8.11 presents (a1, . . . , a4)obtained by the dynamical system
(8.2.2). We start without actuation and switched on the actuation after
10 seconds. First, we see a1and a2reaching the natural attractor. After
switching on the actuation, a1and a2decrease to zero and a3and a4arises
to the actuated attractor.
Figure 8.11: The state a(a1: top left, a2: top right, a3: bottom left and a4:
bottom right) gained by the dynamical system (8.2.2), starting without an
actuation. After 10 seconds, we switched on the actuation.
This figure shows that the dynamical system 8.2.2 represents the behavior
of the natural and the actuated flow as expected.
8.3 Computation of lift
Based on the dynamical system ( ˙ai)i(8.2.2), we calculate the lift by the
following ansatz with unknown coefficients cij
CL(a1, a2, a3, a4) = c00 +
4
X
i=1
c1iai+c20(An)2+c40(An)4.(8.3.1)
There is no limitation on the energy Aaof the actuated case with respect
to increasing B, hence Aais not included in (8.3.1). The ansatz (8.3.1) is
129
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
motivated by a global momentum balance equation and the constant and
linear terms in (8.3.1) are related to this equation. The lift effect of base-
flow variation can be lumped in c20(An)2+c40(An)4assuming slow transients,
see [90]. We obtain the parameters cij by a least-squares fit of CL((an
i)) and
CL((aa
i)) to the original lift values of the URANS simulation. The (an
i)
are the filtered coefficients of the unactuated case and (aa
i)are the filtered
coefficients of the actuated case. Our goal is to fit the parameters in the
sense that the simulated lift values are reproduced by the lift formula (8.3.1)
in the unactuated case CLU based on (an
i)as well as the lift values in the
actuated case CLA based on (aa
i). This leads to the problem
min
cij
F(cij) = kCLA(·)−CL(cij)((aa
i(·)))k2+kCLN (·)−CL(cij)((an
i(·)))k2.
Finally, after the coefficients cij have been determined the optimization prob-
lem is formulated as
max
B>0CL(a1, a2, a3, a4)(8.3.2)
subject to the ODE system (8.0.4).
Figure 8.12 presents the original mean lift by the URANS simulation
(continuous lines) compared with the calculated mean lift by (8.3.1) (dashed
lines).
8.4 Numerical investigation
We consider in this section a high-lift configuration with observation region
Ωpresented in Figure 7.1, see Setting 7.1 and [66] for more details.
The actuation amplitude to determine the set of the actuated snap-
shots was B= 3.5888 and we worked with the parameters σn= 0.5906,
σa=−2.0042,ωn= 5.5407 and ωa= 14.8412. Analogously to (8.2.5), we
calculated g31 = 0.0284,g32 = 0.0000,g41 = 0.0000 and g42 =−0.0019. The
parameters β1= 14.75 and β2= 654.0806 are calculated by (8.2.4) with an
actuation amplitude of B∗= 1.19.
Calibrating the parameters cij of the lift formula (8.3.1) to this data, we
get c10 = 2.2238,c11 = 0.2295,c12 =−0.6858,c13 = 1.6717,c14 =−0.2963,
c20 =−8.3606,c40 = 39.7410. Figure 8.13 shows the agreement of CL(ai)
with the lift-values of the URANS simulation, where the aiare the filtered
coefficients.
The mean values differ in both cases, the unactuated and the actuated
one, between the original lift values and the values calculated with CLnot
more than 0.5%, see Figure 8.12. This is negligible, because in contrast to our
130
8.4. NUMERICAL INVESTIGATION
Figure 8.12: Comparison of the original mean lift values by the URANS
simulation (continuous lines) and the mean lift computed by (8.3.1) (dashed
lines).
Figure 8.13: Comparison of the lift values of the URANS simulation (contin-
uous lines) with those obtained by the lift formula based on the filtered co-
efficients a1,··· , a4(dashed lines): Natural flow (left), actuated flow (right).
stationary case without turbulence, we achieve a lift gain of more than 14% in
the full problem with an actuation in contrast to the case without actuation.
Evaluating CLwith the ai’s as the solutions of the dynamical system, once
computed with B= 0 and once with the full actuation B= 3.5888, we get
131
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
mean values of around 1.96 respectively 2.24 and the results presented in
Figure 8.14.
Figure 8.14: Comparison of the lift values of the URANS simulation (contin-
uous lines) with those obtained by the lift formula based on a1,··· , a4of the
dynamical system (dashed lines): Natural flow (left), actuated flow (right).
Solving this dynamical system with several actuation amplitudes B= 0
to B= 3.6, we resolve the average lift values presented in Figure 8.15; for
B= 0 an average lift of 1.96 and for B= 3.5888 an average lift of 2.22.
The optimization problem (8.3.2) yields a lift gain of more than 13%. The
maximal lift is achieved at an actuation amplitude of around Bopt = 2.4,
agreeing with the results of the URANS simulation.
132
8.4. NUMERICAL INVESTIGATION
Figure 8.15: Comparison of the lift coefficients calculated by CL(continuous
lines) with those obtained by the URANS simulation (dashed lines) with the
Wilcox98 k−ω−turbulence model for some chosen actuation amplitudes B.
133
CHAPTER 8. REDUCED-ORDER MODEL (ROM)
134
Chapter 9
The optimal control problem
In this chapter, we want to investigate an optimal control problem (PN)based
on the optimization problem (8.3.2) of the last chapter. Our goal is to reach
as much lift as possible where the actuation amplitude Bis the optimization
parameter.
Let ˆ
B∈Rdenote an upper bound and 0the natural lower bound for the
actuation amplitude B,a10, a20, a30, a40 the initial values for the state a(t),
T∆=Te−Taand g(t) = Bcos(ωat). Then the optimization problem (PN)
looks with a(t) = (a1(t), a2(t), a3(t), a4(t)) as follows:
min JN(a(t), B) := 1
T∆
Te
Z
Ta
−CL(a1(t), a2(t), a3(t), a4(t)) dt+αN
2B2
subject to the reduced order model
˙a1(t) = ˜σna1(t)−ωna2(t)
˙a2(t) = ωna1(1) + ˜σna2(t)
˙a3(t) = ˜σaa3(1) −ωaa4(t) + g31g(t) + g32 ˙g(t)
˙a4(t) = ωaa3(t) + ˜σaa4(t) + g41g(t) + g42 ˙g(t)
a1(0) = a10
a2(0) = a20
a3(0) = a30
a4(0) = a40
(9.0.1)
with the amplifier rates
˜σn(t) = σn−β1(An(t))2−β2(Aa(t))2
˜σa(t) = σa,
135
CHAPTER 9. THE OPTIMAL CONTROL PROBLEM
the control constraint
B∈Bad := {B∈R: 0 ≤B≤ˆ
B}.
Note that it seems as if this system is linear, but it is nonlinear due to ˜σn
with An(t) = pa2
1(t) + a2
2(t)and Aa(t) = pa2
3(t) + a2
4(t).
In the next section, we introduce the first-order optimality system, which
we will need for the numerical investigation.
9.1 First-order necessary optimality conditions
We directly apply the formal Langrange technique to this problem (PN)to
derive the first-order necessary optimality conditions. We do this in a formal
way without considering the exact function spaces providing the background
of this theory.
Following the Lagrange technique, we want to substitute the ODEs in
(9.0.1) by Lagrange multiplication functions λ(t)=(λ1(t), λ2(t), λ3(t), λ4(t))
while the box constraint B∈Bad for the control parameter and the initial
values are not eliminated.
Then, we obtain the Lagrange-functional LN(a(t), B, λ(t)) with the state
variable a(t), the adjoint state λ(t)and the actuation amplitude B:
LN(a, B, λ) := 1
T∆
Te
Z
Ta
−CL(a)dt+αN
2B2−
Te
Z
Ta
(˙a1−˜σna1+ωna2)λ1dt
−1
T∆
Te
Z
Ta
(˙a2−ωna1−˜σna2)λ2dt
−1
T∆
Te
Z
Ta
(˙a3−σaa3+ωaa4−g31g−g32 ˙g)λ3dt
−1
T∆
Te
Z
Ta
(˙a4(t−ωaa3−σaa4−g41g−g42 ˙g)λ4dt.
Following the Lagrange principle, the optimal control ¯
Btogether with
the associating optimal state ¯a(t)has to fulfill the necessary first-order opti-
mality condition of the problem including the minimization of the Lagrange
functional LNwith respect to a(t), B and the box constraint B∈Bad for the
control but without the state equations (9.0.1).
136
9.1. FIRST-ORDER NECESSARY OPTIMALITY CONDITIONS
Thus, the Lagrange-functional LNhas to fulfill
∂LN
∂a (¯a, ¯
B, λ)h= 0 for all possible h(·)with h(0) = 0 (9.1.1)
and with respect to the control parameter Bthe variational inequality
∂LN
∂B (¯a, ¯
B, λ)(B−¯
B)≥0for all B∈Bad.(9.1.2)
The equality (9.1.1) is with h= (h1, h2, h3, h4)equivalent to
∂LN
∂a (¯a, ¯
B, λ)(h) =
Te
Z
Ta
(
4
X
i=1
(c1ihi)+2c20a1h1+ 2c20a2h2+ 4c40(An)2a1h1
+ 4c40(An)2a2h2)dt
−
Te
Z
Ta
(˙
h1−2(β1a1h1+β1a2h2+β2a3h3+β2a4h4)a1
−σnh1+ωnh2−(β1a2
2+β2a2
3+β2a2
4)h1)λ1dt
−
Te
Z
Ta
(˙
h2−2(β1a1h1+β1a2h2+β2a3h3+β2a4h4)a2
−ωnh1−σnh2−(β1a2
1+β2a2
3+β2a2
4)h2)λ2dt
−
Te
Z
Ta
(˙
h3−˜σah3+ωah4)λ3dt
−
Te
Z
Ta
(˙
h4−ωah3−˜σah4)λ4dt
137
CHAPTER 9. THE OPTIMAL CONTROL PROBLEM
and after integration by parts, we derive
∂LN
∂a (¯a, ¯
B, λ)(h) =
Te
Z
Ta
((c11 + 2c20a1+ 4c40(An)2a1)h1+ (c13)h3
+ (c12 + 2c20a2+ 4c40(An)2a2)h2+ (c14)h4)dt
−
Te
Z
Ta
((−˙
λ1−σn−2(β1a1λ1+β1a1λ2)a1
−(β1a2
2+β2a2
3+β2a2
4) + ωnλ2)h1)dt
−
Te
Z
Ta
((−˙
λ2−ωnλ1−σn−2(β1a2λ1+β1a2h2)a2
−(β1a2
1+β2a2
3+β2a2
4))h2)dt
−
Te
Z
Ta
((−˙
λ3+β2a3(λ1+λ2)−˜σaλ3+ωaλ4)h3)dt
−
Te
Z
Ta
((−˙
λ4+β2a4(λ1+λ2)−ωaλ3−˜σaλ4)h4)dt
−
Te
Z
Ta
(h1(T)λ1(T) + h2(T)λ2(T) + h3(T)λ3(T)
+h4(T)λ4(T))dt.
Since h(T)and h(·)can be arbitrarily, λ= (λ1, λ2, λ3, λ4)is the weak solution
of
−˙
λ1−σn−2(β1a1λ1+β1a1λ2)a1+ωnλ2−(β1a2
2+β2a2
3+β2a2
4) = f1
−˙
λ2−ωnλ1−σn−2(β1a2λ1+β1a2λ2)a2−(β1a2
1+β2a2
3+β2a2
4) = f2
−˙
λ3+β2a3(λ1+λ2)−˜σaλ3+ωaλ4=f3
−˙
λ4+β2a4(λ1+λ2)−ωaλ3−˜σaλ4=f4
λ1(T) = λ2(T) = λ3(T) = λ4(T) = 0
(9.1.3)
with f1=c11 + 2c20a1+ 4c40(An)2a1
f2=c12 + 2c20a2+ 4c40(An)2a2
f3=c13
f4=c14
138
9.2. NUMERICAL INVESTIGATION
which we define as the adjoint system. The solution is the adjoint state λ.
The second requirement (9.1.2) leads to the variational inequality
∂LN
∂B (B−¯
B) = αN¯
B(B−¯
B) + 1
T∆
Te
Z
Ta
[(g31 cos(ωat)−g32ωasin(ωat))λ3
+ (g41 cos(ωat)−g42ωasin(ωat)λ4)](B−¯
B)dt≥0
for all ˆ
B≥B≥0. The pointwise analysis of this inequality leads to the
standard projection formula
¯
B=P[0,ˆ
B]{1
αNT∆
Te
Z
Ta
(g31 cos(ωat)−g32ωasin(ωat))λ3
+ (g41 cos(ωat)−g42ωasin(ωat)λ4)dt}.
(9.1.4)
9.2 Numerical investigation
Now, let us research the optimization problem (PN)numerically based on the
optimality system. In this case, we decided to use the gradient-projection
method, see 5.2.2, because, as mentioned before, COMSOL and the integral
term (9.1.4) don’t fit together. One can see at the end of this chapter that
we need about 5 iterations to get the optimal actuation amplitude. That
means that we have to solve both the nonlinear state equations (9.0.1) and
the linear adjoint system (9.1.3) 5 times. To approximate the optimal the
optimal actuation amplitude ¯
Bby just solving the state equation with dif-
ferent amplitudes would take probably more iterations of the nonlinear state
equation, due to the fact that the upper bound ˆ
Bis free to select. For our
optimization problem (PN)the algorithm reads as follows with given Bn:
S1 Calculate an= ((a1)n,(a2)n,(a3)n,(a4)n)as the solution of (9.0.1)
with the current control Bn.
S2 Calculate the adjoint λn= ((λ1)n,(λ2)n,(λ3)n,(λ4)n)from (9.1.3)
with the current state an.
S3 The updated descent direction is
Dn=αNBn+
T
Z0
(g31 cos(ωat)−g32ωasin(ωat))λ3
+ (g41 cos(ωat)−g42ωasin(ωat)λ4)dt∈R.
139
CHAPTER 9. THE OPTIMAL CONTROL PROBLEM
S4 Calculate the stepsize snfrom
min
s>0f(P[0,ˆ
B]{Bn+sDn}).
S5 The updated control Bn+1 is
Bn+1 := P[0,ˆ
B]{Bn+snDn}.
Set n:=n+1and goto S1.
Let us now investigate the Setting 7.1 on page 90 with the parameters
σn, σa, ωn, ωa, g31, g32, g41, g42, β1, β2, c10, c11, c12, c13, , c14, c20
and c40 as selected Section in 8.4 numerically. Therefore, we chose αN=
0.1and we decide to optimize this problem in the time interval [Ta, Te] =
[56.8612,75.8150].
The first reason of this interval is that we want to optimize the lift and
not the transient oscillation, so we select Ta6= 0. Additionally, we have to
choose Taand Tein the way that they are multiples of both wavelengths, the
natural and the actuated one.
We started the gradient-projection method with initial values
((a1)0,(a2)0,(a3)0,(a4)0)=(rn,0,0,0), where rndenotes the radius of the
natural attractor, B0= 0.5and a mesh size of 0.0132 in the time direction.
The optimal calculated actuation amplitude is Bopt = 2.2573 with an
associated averaged lift coefficient of
1
T∆
Te
Z
Ta
CL((a1)opt,(a2)opt,(a3)opt,(a4)opt)dt= 2.2238,
see Figure 9.1, and JN(aopt, Bopt) = −1.9690 as the value of the cost func-
tional. The calculated optimal actuation amplitude Bopt differs slightly from
the optimal value in Figure 8.15, due to the term αN
2B2in the cost func-
tional. The Figure 9.1 presents the optimal lift coefficient over the time
interval [Ta, Te]and the Figures 9.2 and 9.3 the optimal states (a1,··· , a4)
and the associating adjoint (λ1,··· , λ4), respectively.
Unfortunately, we have no simulations of the full k−ωturbulence system
with our optimal actuation amplitude Bopt.The simulation with an amplitude
nearest to Bopt we have is a simulation with an actuation amplitude of 2.39.
The results of the lift values are compared in Figure 9.5.
140
9.2. NUMERICAL INVESTIGATION
Figure 9.1: The lift coefficient CLof the calculated optimal state over the
interval [Ta, Te].
Figure 9.2: The calculated optimal state (a1: top left, a2: top right, a3:
bottom left, a4: bottom right).
141
CHAPTER 9. THE OPTIMAL CONTROL PROBLEM
Figure 9.3: The calculated adjoint state (λ1: top left, λ2: top right, λ3:
bottom left, λ4: bottom right).
Figure 9.4: The cost functional JNover the number of iterations.
142
9.2. NUMERICAL INVESTIGATION
Figure 9.5: The cost functional JNover the number of iterations.
143
CHAPTER 9. THE OPTIMAL CONTROL PROBLEM
144
Chapter 10
Conclusion
In this thesis, we considered two settings of optimal control problems for
high-lift configurations. In the case of steady-state Navier-Stokes equations
with low Reynolds number, we established first-order necessary optimality
conditions for a problem with an integral state constraint on the drag. The
main theoretical difficulty was the appearance of low regularity controls in
a Dirichlet boundary condition. Afterwards, we considered the second-order
sufficient optimality conditions for the infinite and a finite-dimensional con-
trol space. The optimal control is obtained by direct numerical solution of
the established optimality system and by an SQP-method, where the integral
state constraint was handled by a Penalty term in the cost functional.
An associated nonstationary case with high Reynolds number was inves-
tigated with a WILCOX98 turbulence model. To handle the problem of the
high dimension, a robust reduced order model (ROM) was established fitting
best to the snapshots computed by the full system in the natural and the
actuated state. The ROM reproduces the nonlinear behavior of the system
sufficiently well so that the subsequently optimization problem of periodic
actuation leads to reasonable results. We are now able to solve our optimal
control problem in about 20 minutes by 5 iterations and 4 minutes for one
forward and adjoint system together. Without the model reduction, just one
forward iterations would take about 4 days.
In particular, the application of trust-region proper orthogonal decompo-
sition (TRPOD) could be considered to develop an improved reduced-order
model. In [15] a ROM was used to minimize the total mean drag for a circu-
lar cylinder wake flow by updating the ROM during a (TRPOD) approach,
we refer also to [14].
145
CHAPTER 10. CONCLUSION
146
Chapter 11
Zusammenfassung
In dieser Arbeit haben wir zwei Optimalsteuerungsprobleme für Hochauftrieb-
skonfigurationen untersucht.
Im ersten Fall, der stationären Navier-Stokes Gleichungen mit Kontroll-,
integralen Zustandsbeschränkungen und kleinen Reynolds-Zahlen, haben wir
zunächst die notwendigen Optmalitätsbedingungen erster Ordnung aufgestellt
um das Problem numerisch zu untersuchen. Dabei war die gewünschte
niedrige Regularität der Dirichlet-Randsteuerungen das größte theoretische
Problem. Anschliessend haben wir die hinreichenden Optimalitätsbedingun-
gen zweiter Ordnung für unendlich und endlich dimensionale Steuerungen
aufgestellt. Numerisch haben wir das Problem einerseits als direkte Lösung
des Optimalitätssystems und andererseits mit Hilfe der SQP-Methode un-
tersucht. Zum Abschluss dieses Themenbereichs wurde noch die Konvergenz
der SQP-Methode bewiesen.
Der instationäre Fall wurde mit grossen Reynolds-Zahlen und zugehöri-
gen Turbulenzen betrachtet. Die Turbulenzen wurden durch das WILCOX98
Modell beschrieben, was zu einem riesigen Rechenaufand führt. Alleine eine
Vorwärtsrechnung der Zustandsgleichung hat bei vergleichbaren Problemen
mehr als 4 Tage gedauert. Zur Lösung dieses Problems haben wir eine Mod-
ellreduktion durchgeführt und ein reduced-order model (ROM) aufgestellt,
welches am besten zu vorher berechneten Snapshots passt. Wir haben es
geschafft, dass dieses ROM die nichtlineare Struktur des Systems hinreichend
gut widerspiegelt, so dass eine Optimierung auf dessen Basis möglich ist
und sinnvolle Resultate liefert. Desweiteren gelang es uns das Optimals-
teuerungsproblem innerhalb von etwa 20 Minuten zu lösen, bei 5 Iterationen
und 4 Minuten Dauer für eine Vorwärts- und eine adjungierte Gleichung.
147
CHAPTER 11. ZUSAMMENFASSUNG
148
Bibliography
[1] F. Abergel and E. Casas. Some optimal control problems associated to
the stationary Navier-Stokes equation. In Mathematics, climate and en-
vironment (Madrid, 1991), volume 27 of RMA Res. Notes Appl. Math.,
pages 213–220. Masson, Paris, 1993.
[2] F. Abergel and R. Temam. On some control problems in fluid mechan-
ics. Theoret. Comput. Fluid Dynam., 1:303–325, 1990.
[3] R. A. Adams. Sobolev Spaces. Academic Press, Boston, 1978.
[4] K. Afanasiev and M. Hinze. Adaptive control of a wake flow using
proper orthogonal decomposition. In Lecture Notes in Pure and Applied
Mathematics 216, 317-332, Shape Optimization and Optimal Design,
Marcel Dekker.
[5] W. Alt. The Lagrange-Newton method for infinite-dimensional op-
timization problems. Numer. Funct. Anal. and Optim., 11:201–224,
1990.
[6] W. Alt. Parametric optimization with applications to optimal con-
trol and sequential quadratic programming. Bayreuther Mathematische
Schriften, 1991.
[7] W. Alt. Sequential quadratic programming in Banach spaces. In Ad-
vances in Optimization, volume 382 of Lecture Notes in Economics
and Mathematical systems, pages 281–301, New York, 1992. Springer–
Verlag.
[8] W. Alt. The Lagrange-Newton method for infinite-dimensional opti-
mization problems. Control and Cybernetics, 23:87–106, 1994.
[9] W. Alt. Discretization and mesh independence of Newton’s method for
generalized equations. Preprint, 1997.
149
BIBLIOGRAPHY
[10] W. Alt and K. Malanowski. The Lagrange-Newton Method for Non-
linear Optimal Control Problems. Computational Optimization and
Applications, 2:77–100, 1993.
[11] W. Alt and K. Malanowski. The Lagrange-Newton method for state
constrained optimal control problems. Comp. Optimization and Appl.,
4(3):217–239, 1995.
[12] C. Amrouche and V. Girault. On the existence and regularity of the
solution of Stokes problem in arbitrary dimension. Proc. Japan Acad.
Ser. A Math. Sci., 67(5):171–175, 1991.
[13] Raymond J.-P.. Arada, N. and F. Tröltzsch. On an augmented La-
grangian SQP method for a class of optimal control problems in Ba-
nach spaces. Computational Optimization and Applications, 22:369–
398, 2002.
[14] E. Arian, M. Fahl, and E.W. Sachs. Trust-region proper-orthogonal
decomposition for flow control. Technical Report 2000-25, Universität
Trier, 2000.
[15] M. Bergmann and L. Cordier. Optimal control of the cylinder wake
in the laminar regime by trust-region methods and pod reduced-order
models. J. Comput. Phys., 227(16):7813–7840, 2008.
[16] M. Bergmann, L. Cordier, and J.-P.. Brancher. Optimal rotary control
of the cylinder wake using proper orthogonal decomposition reduced
order model. Phys. Fluids, 17:1–21, 2005.
[17] T. Bewley, P. Moin, and R. Temam. DNS-based predictive control of
turbulence: an optimal benchmark for feedback algorithms. J. Fluid
Mech., 447:179–225, 2001.
[18] M. Braack and T. Richter. Solutions of 3d Navier-Stokes bench-
mark problems with adaptive fnite elements. Computers and Fluids.,
35(4):372–392, 2006.
[19] A. Carnarius, B. Günther, F. Thiele, D. Wachsmuth, F. Tröltzsch,
and J.C. Reyes. Numerical study of the optimization of separation
control. In Proceedings of the 45th AIAA Aerospace Sciences Meeting
and Exhibit, Reno, 8-11 January 2007, AIAA 2007-58.
150
BIBLIOGRAPHY
[20] E. Casas, M. Mateos, and J.-P. Raymond. Error estimates for the nu-
merical approximation of a distributed control problem for the steady-
state Navier-Stokes equations. SIAM J. Control Optim., 46(3):952–982,
2007.
[21] E. Casas and F. Tröltzsch. Second-order necessary and sufficient opti-
mality conditions for optimization problems and applications to control
theory. SIAM J. Optim., 13(2):406–431, 2002.
[22] E. Casas, F. Tröltzsch, and A. Unger. Second order sufficient optimality
conditions for some state-constrained control problems of semilinear
elliptic equations. SIAM J. Control and Optimization, 38:1369–1391,
2000.
[23] Eduardo Casas. The Navier-Stokes equations coupled with the heat
equation: analysis and control. Control Cybernet., 23(4):605–620, 1994.
Modelling, identification, sensitivity analysis and control of structures.
[24] Eduardo Casas. An optimal control problem governed by the evolution
Navier-Stokes equations. In Optimal control of viscous flow, pages 79–
95. SIAM, Philadelphia, PA, 1998.
[25] L. Cattabriga. Su un problema al contorno relativo al sistema di
equazioni di Stokes. Rend. Sem. Mat. Univ. Padova, 31:308–340, 1961.
[26] R. Dautray and J.L. Lions. Mathematical Analysis and Numerical
Methods for Science and Technology. Springer-Verlag, address = Berlin,
year = 2000.
[27] R. Dautray and J.L. Lions. Mathematical Analysis and Numerical
Methods for Science and Technology, Vol. 6. Springer Verlag, Berlin,
2000.
[28] J. C. de los Reyes and K. Kunisch. A semi-smooth Newton method
for control constrained boundary optimal control of the Navier-Stokes
equations. Nonlinear Anal., 62(7):1289–1316, 2005.
[29] J. C. de los Reyes and F. Tröltzsch. Flow control with regularized
state constraints. in: Active Flow Control, Notes on Numerical Fluid
Mechanics and Multidisciplinary Design (NNFM), 95:353–366, 2007.
[30] J. C. de los Reyes and I. Yousept. Regularized state-constrained bound-
ary optimal control of the navier-stokes equations. Journal of Mathe-
matical Analysis and Applications., 356:257–279, 2009.
151
BIBLIOGRAPHY
[31] J.C. de los Reyes, P. Merino, J. Rehberg, and F. Tröltzsch. Optimal-
ity conditions for state-constrained pde control problems with time-
dependent controls. To appear in Control and Cybernetics.
[32] P. Deuflhard. Newton Methods for Nonlinear Problems. Affine Invari-
ance and Adaptive Algorithms. Band 35 der Reihe Springer Series in
Computational Mathematics. Springer-Verlag, Berlin, 2004.
[33] A. L. Dontchev. Uniform convergence of the Newton method for Aubin
continuous maps. Serdica Math. J., 22(3):385–398, 1996.
[34] A. L. Dontchev, W. W. Hager, A. B. Poore, and B. Yang. Optimality,
stability, and convergence in nonlinear control. Applied Math. and
Optimization, 31:297–326, 1995.
[35] Asen L. Dontchev. Implicit function theorems for generalized equa-
tions. Math. Program., 70(1):91–106, 1995.
[36] Asen L. Dontchev. Local analysis of a Newton-type method based on
partial linearization. volume 32, pages 295–306, 1996.
[37] Asen L. Dontchev and William W. Hager. Implicit functions, lipschitz
maps, and stability in optimization. Math. Oper. Res., 19(3):753–768,
1994.
[38] J. Dušek, P. Le Gal, and P. Fraunié. A numerical und theoretical study
of the Hopf bifurcation in a cylinder wake. J. Fluid Mech., 264:59–80,
1994.
[39] R. Farwig, G. P. Galdi, and H. Sohr. A new class of weak solutions
of the Navier-Stokes equations with nonhomogeneous data. J. Math.
Fluid Mech., 8(3):423–444, 2006.
[40] R. Farwig, G. P. Galdi, and H. Sohr. Very weak solutions and large
uniqueness classes of stationary Navier-Stokes equations in bounded
domains of R2.J. Differential Equations, 227(2):564–580, 2006.
[41] C.A.J. Fletcher. Computational Galerkin Methods. Springer-Verlag,
1984.
[42] R. Fletcher. Practical Methods of Optimization. Wiley, 1987.
[43] M. Gad-el Hak. Flow Control. Passive, Active, and Reactive Flow
Management. Cambridge University Press, 2000.
152
BIBLIOGRAPHY
[44] G. P. Galdi. An introduction to the mathematical theory of the Navier-
Stokes equations. Springer, New York, 1992.
[45] G. P. Galdi, C. G. Simader, and H. Sohr. A class of solutions to
stationary Stokes and Navier-Stokes equations with boundary data in
W−1/q,q.Math. Ann., 331(1):41–74, 2005.
[46] P. E. Gill, W. Murray, and M. H. Wright. Practical Optimization.
Academic Press, London, 1981.
[47] V. Girault and P. A. Raviart. Finite Element Methods for Navier Stokes
Equations. Springer Verlag Berlin, 1986.
[48] H. Goldberg and F. Tröltzsch. On a Lagrange-Newton Method for a
nonlinear parabolic boundary control problem. Optimization Methods
and Software, 8:225–247, 1998.
[49] W.R. Graham, J. Peraire, and K.Y. Tang. Optimal control of vortex
shedding using low-order models. Part I:Open-loop model development.
Int. Num. Meth. Eng., 44:945–972, 1999.
[50] R. Griesse and J. C. de los Reyes. State-constrained optimal control
of the three-dimensional stationary Navier-Stokes equations. J. Math.
Anal. Appl., 343(1):257–272, 2008.
[51] M. D. Gunzburger, L. S. Hou, and Th. P. Svobodny. Analysis and finite
element approximation of optimal control problems for the stationary
Navier-Stokes equations with Dirichlet controls. RAIRO Modél. Math.
Anal. Numér., 25(6):711–748, 1991.
[52] M. D. Gunzburger and S. Manservisi. The velocity tracking problem for
Navier-Stokes flows with boundary control. SIAM J. Control Optim.,
39:594–634, 2000.
[53] M. Hintermüller and M. Hinze. A SQP-Semi-Smooth Newton-type
Algorithm applied to Control of the instationary Navier-Stokes System
Subject to Control Constraints. Siam J. Optim., 16:1177–1200, 2006.
[54] M. Hinze and K. Kunisch. Second-order methods for optimal control of
time-dependent fluid flow. SIAM J. Control Optim., 40:925–946, 2001.
[55] M. Hinze and K. Kunisch. Second order methods for boundary control
of the instationary Navier-Stokes system. ZAMM, 84:171–187, 2004.
153
BIBLIOGRAPHY
[56] P. Holmes, J. L Lumley, and G. Berkooz. Turbulence, Coherent Struc-
tures, Dynamical Systems and Symmetry. Cambridge University Press,
Cambridge, 1998.
[57] C. John, B. R. Noack, M. Schlegel, F. Tröltzsch, and D. Wachsmuth.
Optimal Boundary Control Problems Related to High-Lift Configura-
tions. Active Flow Control II. Notes on Numerical Fluid Mechanics
and Multidisciplinary Design, 108/2010.
[58] C. John and D. Wachsmuth. Optimal Dirichlet boundary control of
Navier-Stokes equations with state constraint. Numerical Functional
Analysis and Optimization, 30(11&12).
[59] B.H. Jørgensen, J.N. Sørensen, and M. Brøns. Low-dimensional mod-
elling of a driven cavity flow with two free parameters. Theoret. Com-
put. Fluid Dynamics, 16:299–317, 2003.
[60] N. H. Josephy. Newton’s method for generalized equations. Technical
report, 1979.
[61] K. Kunisch and B. Vexler. Constrained Dirichlet boundary control
in L2for a class of evolution equations. SIAM J. Control Optim.,
46(5):1726–1753, 2007.
[62] K. Kunisch and S. Volkwein. Galerkin proper orthogonal decomposition
methods for a general equation in fluid dynamics. SIAM Journal on
Numerical Analysis, 40:492–515, 2002.
[63] K. Kunisch and L. Xie. Suboptimal feedback control of
flow separation by POD model reduction. SIAM, DOI:
10.1137/1.9780898718935.ch12, 2006.
[64] F.-S. Kupfer and E. W. Sachs. Numerical solution of a nonlinear
parabolic control problem by a reduced SQP method. Computational
Optimization and Applications, 1:113–135, 1992.
[65] O.A. Ladyzhenskaya. The Mathematical Theory of Viscous Incompress-
ible Flow. Gordan and Breach, 1963.
[66] D.M. Luchtenburg, B. Günther, B.R. Noack, R. King, and G. Tadmor.
A generalized mean-field model of the natural and high-frequency actu-
ated flow around a high-lift configuration. J. Fluid Mech., 623:283–316,
2009.
154
BIBLIOGRAPHY
[67] K. Malanowski. Two-norm approach in stability and sensitivity analysis
of optimization and optimal control problems. Adv. Math. Sci. Appl.,
2(2):397–443, 1993.
[68] K. Malanowski. Sufficient optimality conditions in optimal control.
Technical report, 1994.
[69] E. Marušić-Paloka. Solvability of the Navier-Stokes system with L2
boundary data. Appl. Math. Optim., 41(3):365–375, 2000.
[70] H. Maurer and J. Zowe. First- and second-order conditions in infinite-
dimensional programming problems. Math. Programming, 16:98–110,
1979.
[71] V. Maz’ya and J. Rossmann. Mixed boundary value problems for the
Navier-Stokes system in polyhedral domains. Mathematical Physics,
arXiv:math-ph/0602054v1, 2006.
[72] V. Maz’ya and J. Rossmann. Lpestimates of solutions to mixed bound-
ary value problems for the Stokes system in polyhedral domains. Math.
Nachr., 280(7):751–793, 2007.
[73] I. Neitzel, U. Prüfert, and T. Slawig. Strategies for time-dependent pde
control with inequality constraints using an integrated modeling and
simulation environment. Numericl Algorithms, 50(3):241–269, 2009.
[74] B.R. Noack, K. Afanasiev, M. Morzyński, G Tadmor, and F. Thiele.
A hierarchy of low-dimensional models for the transient and post-
transient cylinder wake. J. Fluid Mech., 497:335–363, 2008.
[75] B.R. Noack and G.S. Copeland. On a stability property of ensemble-
averaged flow. Tech. report 03/2000. Hermann-Föttinger-Institut für
Strömungsmechanik., 3, 2000.
[76] B.R. Noack, P. Papas, and P.A. Monkewitz. The need for a pressure-
term representation in empirical Galerkin models of incompressible
shear flows. J. Fluid Mech., 523:339–365, 2005.
[77] B.R. Noack, G. Tadmor, and M. Morzyński. Actuation models and
dissipative control in empirical Galerkin models of fluid flows. In In
Proceedings of the 2004 American Control Conference, pages 5722–
5727, Daytona, OH, USA, 2004. American Automatic Control Council
(AACC).
155
BIBLIOGRAPHY
[78] J. Nocedal and S. J. Wright. Numerical Optimization. Springer-Verlag,
New York, 1999.
[79] J.-P. Raymond. Stokes and Navier-Stokes equations with nonhomoge-
neous boundary conditions. Ann. Inst. H. Poincaré Anal. Non Linéaire,
24(6):921–951, 2007.
[80] J. C. de los Reyes and K. Kunisch. A semi-smooth Newton method
for regularized state-constrained optimal control of the Navier-Stokes
equations. Computing, 78(4):287–309, 2006.
[81] J. C. de los Reyes and F. Tröltzsch. Optimal control of the stationary
Navier-Stokes equations with mixed control-state constraints. SIAM
J. Control Optim., 46(2):604–629, 2007.
[82] S. M. Robinson. Stability theory for systems of inequalities, part ii:
Differentiable nonlinear systems. SIAM J. Numer. Anal., 1979.
[83] S. M. Robinson. Strongly regular generalized equations. Mathematics
of Operation Research, 1980.
[84] T. Roubíček and F. Tröltzsch. Lipschitz stability of optimal controls
for the steady-state Navier-Stokes equations. Control and Cybernetics,
32(3):683–705, 2002.
[85] C.W. Rowley and V. Juttijudata. Model-based control and estimation
of cavity flow oscillations. In In Proceedings of the 44th IEEE Confer-
ence on Decision and Control and European Control Conference, pages
512–517, Saint-Martin d’Hères, France., 2005. European Union Control
Association (EUCA).
[86] M. Samimy, M. Debiasi, E. Caraballo, A. Serrani, X. Yuan, J. Little,
and J. Myatt. Feedback control of subsonic cavity flows using reduced-
order models. J. Fluid Mech., 579:315–346, 2007.
[87] Günther B. Schatz, M. and F. Thiele. Computational investigation
of separation control for high-lift airfoil flows. Active Flow Control,
Springer, 2006.
[88] M. Schatz and F. Thiele. Numerical Study of High-Lift Flow with
Separation Control by Periodic Excitation. AIAA paper, 2001-0296,
2001.
156
BIBLIOGRAPHY
[89] M. Schatz, F. Thiele, R. Petz, and W. Nitsche. Separation control
by periodic excitation and its application to a high-lift configuration.
AIAA, (2507), 2004.
[90] Luchtenburg D.M. Noack B.R. Aleksic K. Pastoor M. King R. Schlegel,
M. and G Tadmor. Turbulence Control Based on Reduced-Order Mod-
els and Nonlinear Control Design. In Active Flow Control II.
[91] Zh. W. Shen. A note on the Dirichlet problem for the Stokes system
in Lipschitz domains. Proc. Amer. Math. Soc., 123(3):801–811, 1995.
[92] S. Siegel, K. Cohen, and T. McLaughlin. Feedback control of a circular
cylinder wake in experiment and simulation. AIAA Paper, pages 2003–
3569, 2003.
[93] S.G. Siegel, J. Seidel, C. Fagley, D.M. Luchtenburg, K. Cohen, and
T. McLaughlin. Low-dimensional modelling of a transient cylinder
wake using double proper orthogonal decomposition. J. Fluid Mech.,
610:1–42, 2008.
[94] J.T. Stuart. On the non-linear mechanics of hydrodynamics stability.
J. Fluid Mech., 4:1–21, 1958.
[95] J.T. Stuart. Nonlinear stability theory. Annu. Rev.Fluid Mech., 3:347–
370, 1971.
[96] G. Tadmor, O. Lehmann, B. R. Noack, L. Cordier, J. Delville, j. P. Bon-
net, and M. Morzyński. Reduced order models for closed-loop wake
control. Phil. Trans. R. Soc. A, 369:1513 – 1624, 2011.
[97] R. Temam. Navier-Stokes equations. North Holland, Amsterdam, 1979.
[98] F.H. Tinapp and W. Nitsche. On active control of high-lift flow. In In
W. Rodi and D. Laurence, editors, Proc. 4th Int. Symposium on En-
gineering Turbulence Modelling and Measurements, Corsica, Elsevier
Science, 1999.
[99] F. Tröltzsch. Optimality conditions for parabolic control problems and
applications, volume 62 of Teubner Texte zur Mathematik. B.G. Teub-
ner Verlagsgesellschaft, Leipzig, 1984.
[100] F. Tröltzsch. An SQP method for the optimal control of a nonlinear
heat equation. Control and Cybernetics, 23(1/2):267–288, 1994.
157
BIBLIOGRAPHY
[101] F. Tröltzsch. Convergence of an SQP–Method for a class of nonlinear
parabolic boundary control problems. In W. Desch, V. Kappel, and
K. Kunisch, editors, Control and Estimation of Distributed Parame-
ter Systems: Nonlinear Phenomena. Int. Series of Num. Mathematics,
volume 118, pages 343–358, 1994.
[102] F. Tröltzsch. Optimal Control of Partial Differential Equations. The-
ory, Methods and Applications. Graduate Studies in Mathematics, Vol-
ume 112. American Mathematical Society, Providence, 2010.
[103] F. Tröltzsch and D. Wachsmuth. Second-order sufficient optimality
conditions for the optimal control of Navier-Stokes equations. ESAIM:
COCV, 12:93–119, 2006.
[104] A. Unger. Hinreichende Optimalitätsbedingungen 2. Ordnung und Kon-
vergenz des SQP-Verfahrens für semilineare elliptische Randsteuer-
probleme. PhD thesis, Technische Universität Chemnitz, 1997.
[105] S. Volkwein. Proper orthogonal decomposition and singular value de-
composition. SFB-Preprint, 153, 1999.
[106] S. Volkwein. Optimal control of a phase-field model using proper or-
thogonal decomposition. ZAMM, 81:83–97, 2001.
[107] S. Volkwein. Lagrange-SQP techniques for the control constrained op-
timal boundary control problems for the Burgers equation. Computa-
tional Optimization and Applications, 26:253–284, 2003.
[108] S. Volkwein and F. Tröltzsch. Pod a-posteriori error estimates for
linear-quadratic optimal control problems. Computational Optimiza-
tion and Applications, 44:83–115, 2009.
[109] D. Wachsmuth. Optimal control of the unsteady Navier-Stokes equa-
tions. PhD thesis, Technische Universität Berlin, 2006.
[110] D. Wachsmuth. Sufficient second-order optimality conditions for con-
vex control constraints. J. Math. Anal. App., 319:228–247, 2006.
[111] Y. Wang, F. Tröltzsch, and G. Bärwolff. The POD Dirichlet Boundary
Control of the Navier-Stokes Equations: A Low-dimensional Approach
to Optimal Control with High Smoothness. Technical Report 23-2009,
Technische Universität Berlin, 2009.
[112] E. Wassen and F. Thiele. Active control of a model vehicle wake.
Journal of Turbulence, 2006.
158
BIBLIOGRAPHY
[113] D.C. Wilcox. Simulation of Transition with a Two-Equation Turbu-
lence Model. AIAA JOURNAL, 32, 1994.
[114] D.C. Wilcox. Turbulence Modeling for CFD. DCW Industries, La
Cañada, 2006.
[115] J. Zowe and S. Kurcyusz. Regularity and stability for the mathematical
programming problem in Banach spaces. Appl. Math. Optim., 5(1):49–
62, 1979.
159
