Shape optimization of a breakwater
Moritz Keuthen
a,∗
, Daniel Kraft
b
a
Chair of Mathematical Optimization, Department of Mathematics, Technische Universit¨at M¨unchen, Boltzmannstr. 3, 85748 Garching b.
M¨unchen, Germany
b
Institute of Mathematics, University of Graz, NAWI Graz, Universit¨atsplatz 3, 8010 Graz, Austria
Abstract
In this paper we optimize the shape of a breakwater which protects a harbor basin from incoming waves. More specif-
ically, our objective is reducing the harbor resonance due to long range ocean waves. We will consider the complex
valued Helmholtz equation as our model state equation and minimize the average wave height in the harbor basin
with the shape of the breakwater as optimization variable. The geometry will be described by the level set method,
i.e. the domain is given as the sub-zero level set of a function. In contrast to many publications we use the volume
expression of the shape derivative, which lends itself naturally to a level set update via a transport equation. The vol-
ume expression requires less regular finite element functions than the Hadamard form. The model problem naturally
allows to treat geometric constraints in the form of forbidden regions and we will compare different possibilities to
handle them.
Keywords: Shape Optimization, Level Set Method, Helmholtz Equation, Harbor Resonance, Geometric Constraints
1. Introduction
Since the introduction of the level set method to shape optimization at the turn of the century it has developed
into one of the most powerful techniques in this area. Some early works are [1, 2, 3, 4, 5]. Many publications deal
with structural optimization, usually the Ersatz material approach in the hold-all domain D is used to compute the
mechanical properties of the structure and the domain Ω is never explicitly resolved. But there are also approaches
where the domain is exactly meshed, consult for instance [6, 7, 8, 9, 10, 11, 12]. The literature is extensive, we refer
to the review paper [13] for an overview of level set based methods in structural topology and shape optimization.
The survey [14] presents the level set method in combination with inverse problems and optimal design. For a more
comprehensive exposition of the broad field of shape optimization we mention the monographs [15, 16, 17, 18, 19,
20, 21, 22, 23].
In this paper our objective is to reduce the resonance of a harbor due to long range ocean waves. In practice many
publications (e.g. [24, 25, 26, 27, 28]) employ the mild slope equation to model wave effects in coastal areas. It was
first derived by Berkhoff [29] and can be written as
∇ · CC
g
∇ϕ + k
2
CC
g
ϕ = 0, (1)
where ϕ is the horizontal variation in velocity potential, k is the wave number, ω is the wave frequency, C = ω/k
is the wave celerity, C
g
=
C
2
1 +
2kh
sinh 2kh
is the group velocity and h is the water depth. Usually it is enriched with
additional effects such as partial absorption boundaries, bottom friction, entrance loss, etc. In this note we will assume
as simplification that the water depth is constant throughout which leads to the well-known Helmholtz equation
∆ϕ + k
2
ϕ = 0. (2)
∗
Corresponding author
Email addresses: keuthen@ma.tum.de (Moritz Keuthen), daniel.kraft@uni-graz.at (Daniel Kraft)
Preprint submitted to Elsevier November 7, 2014
Note that the potential ϕ in this case is complex valued. To the best of our knowledge such a shape optimization
problem was only briefly considered in the thesis [30]. The author used the real-valued Helmholtz equation as state
equation, an explicit discrete geometry description via the finite element mesh and computed the shape derivative by
differentiating with respect to nodal coordinates.
We will consider the complex valued Helmholtz equation as our model state equation and want to minimize the
average wave height in the harbor basin with the shape of the breakwater as optimization variable. The geometry
will be described by the level set method [31, 32, 33, 34], i.e. the domain Ω is given as the sub-zero level set of a
function Φ : D → R, where D ⊂ R
2
is the hold-all domain. In the numerical realization the level set function is
given on a regular grid and we approximate ∂Ω with one or more polygonal curves. In each iteration we construct
a triangulation of D which resolves the interface and use it to compute the state, adjoint state and gradient of the
cost functional. The model state equation is given on an exterior domain. Since we want to use finite elements we
follow the presentation in [35] and reformulate the problem on a bounded domain. In contrast to many publications
we use the volume expression of the shape derivative, which lends itself naturally to a level set update via a transport
equation. It requires less regular finite element functions and in our experience the volume expression is numerically
more stable than the Hadamard form. This assessment is shared in the recent papers [36, 37]. In [37] the volume
expression of the shape derivative and the level set method are also used.
The model problem naturally involves geometric constraints in the form of forbidden regions. We strictly enforce
those constraints by projecting the gradient onto a suitable admissible set. From a theoretic point of view the scalar
product used for the projection should be at least as smooth as the scalar product used to determine the gradient. In
the numerical experiments we will study different choices of gradient and projection and see that sometimes a less
regular projection leads to better results.
In the following we will use the shorthand notation X for the product space X
d
.
2. Description of the physical model
We study the situation depicted in Figure 1. We have an isle bounded by the contour Γ
L
, some breakwaters given
by Γ
B
and a surrounding ocean denoted by Ω
+
. We want to compute the scattered wave u induced by an incoming
planar monochromatic wave z(x) = e
ikd
T
x
with incident direction d ∈ R
2
and wave number k > 0. The total surface
perturbation is given by y = u + z. Our mathematical model is given by
∆y + k
2
y = 0 in Ω
+
ay +
∂y
∂n
= 0 on Γ
I
:= Γ
B
∪ Γ
L
.
(3)
Here a : D → C describes the absorption coefficient at the boundary. The boundaries Γ
L
, Γ
B
are assumed to be
Lipschitz. Let us discuss appropriate boundary conditions. It is a standard assumption that the scattered wave satisfies
the Sommerfeld radiation condition
iku −
∂u
∂R
= o(R
−
1
2
), for R → ∞.
If we want to study this problem in weak form on the unbounded domain Ω
+
we need to introduce different
weighted Sobolev spaces for the test and ansatz functions and include the Sommerfeld radiation condition in the
ansatz space. See [38, chapter 4] and [35, section 2.3] for more details. This approach however leads to various
difficulties in the numerical realization.
An alternative is to decompose the domain Ω
+
disjointly into a bounded domain Ω and an unbounded domain Ω
a
by introducing an artificial smooth boundary Γ
a
, such that Ω
+
= Ω ∪ Γ
a
∪ Ω
a
. The problem (3) is then equivalent to
the following coupled problem (c.f. [39])
∆y
−
+ k
2
y
−
= 0 in Ω
ay
−
+
∂y
−
∂n
= 0 on Γ
I
y
−
= y
+
on Γ
a
∂y
−
∂n
=
∂y
+
∂n
on Γ
a
∆y
+
+ k
2
y
+
= 0 in Ω
a
iku
+
−
∂u
+
∂R
= o(R
−
1
2
), for R → ∞.
(4)
2
Figure 1: The domain Ω
For a given u
−
on Γ
a
one can solve the unbounded Dirichlet problem (recall u = y −z)
∆u
+
+ k
2
u
+
= 0 in Ω
a
u
+
= u
−
on Γ
a
iku
+
−
∂u
+
∂R
= o(R
−
1
2
), for R → ∞,
compare [38]. If we have the solution u
+
we can easily compute
∂u
+
∂n
=
∂u
−
∂n
on Γ
a
. We denote the mapping u
−
7→
∂u
−
∂n
by G
e
and observe that G
e
∈ L(H
1
2
(Γ
a
), H
−
1
2
(Γ
a
)). This operator is called the Dirichlet-to-Neumann (DtN) operator.
There exists an integral and a series representation of the non-local operator G
e
, the integral variant is also called the
Steklov-Poincar
´
e operator. The Neumann condition
∂y
−
∂n
=
∂y
+
∂n
on Γ
a
can be reformulated by
∂y
−
∂n
=
∂u
−
∂n
+
∂z
∂n
= G
e
u
−
+
∂z
∂n
= G
e
y
−
−G
e
z +
∂z
∂n
.
Now we replace G
e
by some yet unspecified G ∈ L(H
1
2
(Γ
a
), H
−
1
2
(Γ
a
)), which might be some approximation of the
exact solution operator G
e
. Hence we arrive at the bounded problem
∆y + k
2
y = 0 in Ω
ay +
∂y
∂n
= 0 on Γ
I
∂y
∂n
= Gy − Gz +
∂z
∂n
on Γ
a
,
(5)
which is equivalent to (3) for the choice G = G
e
.
Let us fix some conventions and notation at this point. We define the usual bilinear L
2
-scalar product for real-
valued functions on some set A ⊂ R
d
as (·, ·)
L
2
(A)
and the corresponding sesquilinear form as ( f, g)
L
2
C
(A)
:= ( f, g)
L
2
(A)
for some complex valued functions f, g. Furthermore we introduce the real-valued scalar product
( f, g)
L
2
R
(A)
:= Re( f, g)
L
2
C
(A)
.
The norm k·k
L
2
R
(A)
induced by this scalar product coincides with the norm induced by the sesquilinear form. Hence the
elements of the space
L
2
R
(A) := {f : A 7→ C | kf k
L
2
R
(A)
< ∞}
coincide with the elements of {f : A 7→ C | kf k
L
2
C
(A)
< ∞}, but since we use the (·, ·)
L
2
R
(A)
scalar product we have a
different Hilbert space structure. Other Hilbert spaces will be treated analogously (e.g. H
1
R
(A)).
Let us get back to the model problem. We make the following simplifying assumption:
Assumption 1. We choose G as the 0th-order approximation of G
e
hGy, ϕi
H
−
1
2
R
(Γ
a
),H
1
2
R
(Γ
a
)
= (iky, ϕ)
L
2
R
(Γ
a
)
,
(c.f. [35, section 3.]). Furthermore the absorption coefficient is set to a ≡ 0 (i.e. perfect reflection).
3
Remark 1. The choice of G is a rather crude simplification, but in this paper our focus is more on methodology
and not so much on realistic modeling. The low-order approximation of G
e
introduces artificial reflections at Γ
a
, but
since this makes the shape optimization problem presumably harder we accept that for the moment. For more evolved
methods of treating the artificial boundary Γ
a
and the operator G
e
we refer to [35, chapter 3].
Setting the absorption coefficient to zero implies that there is no damping effect by absorption of energy at the
reflecting boundary. We note again that this presumably makes the optimization problem harder since small design
changes might have large non-local effects because of wave interference. Furthermore, while breakwaters are efficient
in reducing the wave amplitude for waves with short wave periods (like wind waves), this does not hold for long
wave periods. Harbor oscillations and resonance due to long waves has been widely studied in coastal engineering
literature, see for example [40] and the references therein.
The weak formulation of (5), with G, a chosen to satisfy Assumption 1, reads
(
Find y ∈ H
1
R
(Ω) :
b(y, ϕ) = f (ϕ), ∀ϕ ∈ H
1
R
(Ω),
(6)
where we define
b(y, ϕ) := (∇y, ∇ϕ)
L
2
R
(Ω)
− k
2
(y, ϕ)
L
2
R
(Ω)
− (iky, ϕ)
L
2
R
(Γ
a
)
,
f (ϕ) := (
∂z
∂n
− ikz, ϕ)
L
2
R
(Γ
a
)
.
(7)
Results concerning the existence of a unique solution of (6) and its regularity are well-known:
Theorem 1. Let Ω be a Lipschitz domain. Then there exists a unique solution y ∈ H
1
R
(Ω) of (6) for any right-hand
side f ∈ H
1
R
(Ω)
∗
and we have kyk
H
1
R
(Ω)
≤ ckf k
H
1
R
(Ω)
∗
for some c > 0.
Proof. We have the Gelfand-triple H
1
R
(Ω) → L
2
R
(Ω) → H
1
R
(Ω)
∗
. Further b(·, ·) is H
1
R
(Ω)-coercive. Hence the
Fredholm alternative holds: Either there exists a unique solution of (6) for any f ∈ H
1
R
(Ω) or there exists a nontrivial
solution y
0
, 0 of the homogenous problem
b(y, ϕ) = 0, ∀ϕ ∈ H
1
R
(Ω).
For our choice of G the solution of the homogenous problem is unique [35, Theorem 3.2].
Theorem 2. Let Ω be a C
2
-domain. The solution y ∈ H
1
R
(Ω) of (6) has the additional regularity y ∈ H
2
R
(Ω).
Proof. Follows from standard regularity results for elliptic equations c.f. [41, Theorem 9.1.20].
3. Shape optimization problem
As announced in the introduction our objective is to minimize the average wave height in the harbor basin Q
(compare Figure 1). Given a solution y of (6) we define the cost functional as
J(Ω, y) =
1
2
kyk
2
L
2
R
(Q)
.
Obviously enclosing the whole harbor basin by a breakwater is not a feasible solution, so we have to introduce an
harbor approach A and demand that Q ∪ A is always part of the ocean. Furthermore we do not want to remove parts
of the island (the inhabitants might complain). Hence we define the admissible set of domains to be
O
ad
:= {Ω ⊂ D | (Q ∪ A) ⊂ Ω, Ω ∩ L = ∅},
where Ω ∈ O
ad
represents the ocean. Of course one could impose additional constraints, e.g. a volume constraint on
the set of admissible domains. We can now formulate the abstract shape optimization problem
minimize J(Ω, y) such that Ω ∈ O
ad
, and y solves (6) on Ω.
In this paper we will not concern ourselves with the question of existence of a solution of the shape optimization
problem. For a detailed discussion of techniques and conditions which guarantee the existence and uniqueness of
solutions to shape optimization problems, we refer to [18, 19]. In order to find a solution (if it exists) we develop
a framework for a gradient descent method. For this we follow the usual optimal control approach of introducing a
design-to-state mapping Ω 7→ y(Ω) and computing the derivative of the reduced objective j(Ω) := J(Ω, y(Ω)). In the
next section we outline the theory of computing shape derivatives and apply it to our problem setting.
4
4. Shape sensitivity analysis
Following the presentation in [42, 18] we introduce a bounded hold-all domain D ⊂ R
d
, D , ∅, which contains
all relevant geometrical objects, and study vector fields V : D → R
d
, satisfying the conditions
V is globally Lipschitz continuous in D, and supp(V) ⊂ D. (8)
Note that in [42, 18] a more general setting is considered, but this suffices for our purposes. If the velocity field V
satisfies condition (8) we associate it with the so called flow map T
t
(V): d/dt
(
T
t
(V)(x)
)
= V(x). It transforms a set
Ω ⊂ D into a new domain
Ω
t
(V) := T
t
(V)(Ω) = {T
t
(V)(x) | ∀x ∈ Ω},
which is also contained in D. If it is clear which velocity field drives the transformation we will drop it to shorten
notation. Furthermore we introduce the spaces of k−times differentiable functions with compact support in D, i.e.
D
k
(D) := {θ ∈ C
k
(D) | supp θ ⊂⊂ D}.
The result [18, Theorem 4.5.1] shows that the transformation T
t
is bi-Lipschitz for some τ > 0 and t ∈ [0, τ] if V
satisfies (8). If V ∈ D
k
(D) then T
t
is a C
k
-diffeomorphism.
Let us recall the classical notions of the Eulerian semiderivative and shape differentiability:
Definition 1. [42, Definition 3.1, Definition 3.2] Let A ⊂ P(D) be an appropriate set of shape variables contained in
the hold-all domain D ⊂ R
d
.
1. The Eulerian semiderivative of j : A → R at Ω ∈ A in a direction V satisfying (8) is given by
d j(Ω; V) := lim
t&0
j(Ω
t
(V)) − j(Ω)
t
,
if the limit exists and is finite.
2. The shape functional j : A → R is said to be shape differentiable at Ω ∈ A, if the Eulerian semiderivative
exists for all V ∈ D
∞
(D) and the map
d j(Ω; ·) : D
∞
(D) → R, V 7→ d j(Ω; V),
is linear and continuous.
Remark 2. 1. The Hadamard-Zolesio structure theorem (c.f. [42, Theorem 3.2]) shows that the support of the
vector distribution d j(Ω; ·) ∈ D
∞
(D)
∗
is contained in ∂Ω ∩ D.
2. The proofs of most formulas for the computation of shape derivatives rely on the representation of the shape
functional j in terms of the transformation T
t
(V). If j depends on the integral over a d-dimensional subset of Ω
t
and a shape depended state y the directional derivative has a natural representation as a volume integral. If the
boundary ∂Ω is smooth enough this can be related via Gauß’ divergence theorem to a boundary representation
in accordance with the structure predicted by the Hadamard-Zolesio theorem.
3. If the boundary ∂Ω of Ω ⊂⊂ D is compact then d j(Ω; ·) is continuous for the D
k
(D)-topology for some k ≥ 0
and d j(Ω; ·) ∈ H
−s
(D) for some s ≥ 0 (c.f. [18, Remark 9.3.1]).
4. The derivative d j(Ω; ·) ∈ D
∞
(D)
∗
is often called the shape gradient. We think that the terms derivative and
gradient should be clearly separated. Whereas the derivative is an element of the dual space of some vector
space, in our terminology the gradient is the Riesz representative of the derivative with respect to some scalar
product. So if d j(Ω; ·) ∈ H
∗
for some Hilbert space H the gradient ∇j(Ω) with respect to the H-scalar-product
is given by
(∇j(Ω), V)
H
= d j(Ω; V), ∀V ∈ H.
5
