# Topology Optimization of Time Dependent Viscous Incompressible Flows

Authors: Mohamed Abdelwahed; Maatoug Hassine
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

In this paper, we consider topology and shape optimization problem related to the nonstationary Navier-Stokes system. The minimization of dissipated energy in the fluid flow domain is discussed. The proposed approach is based on a sensitivity analysis of a design function with respect to the insertion of a small obstacle in the fluid flow domain. Some numerical results show the efficiency and accurate of the proposed approach.


---

## Body

## 1. Introduction

Topological shape optimization in fluid mechanics has wide and valuable applications in hydrodynamic and aerodynamic problems such as the design of car hoods, airplane wings, and inlet shapes for jet engines. Various optimization methods are introduced to determine the optimal design of minimum drag bodies [1–3], diffusers [4], valves [5], and airfoils [6]. The majority of methods dealing with optimal design of flow domains fall into the category of shape optimization and are limited to determine the optimal shape of an existing boundary.It is only recently that topological optimization has been introduced and used in fluid shape optimization problems. It can be used to design features within the domain allowing new boundaries to be introduced into the design. In this context, one of the first approaches is proposed by Borrvall and Petersson in [7]. They implemented the relaxed material distribution approach to minimize the power dissipated in Stokes flow. Later, this approach has been generalized by Guest and Prévost in [8]. They treated the material phase as porous medium where fluid flow is governed by Darcy’s law. In [9], we have proposed a new topological optimization approach, for the Stokes system, based on the topological sensitivity analysis [10–17]. The optimal domain is constructed through the insertion of some obstacles in the initial one. The problem leads to optimizing the obstacles location. The topological sensitivity analysis method consists in studying the variation of the cost function j with respect to the insertion of a small obstacle 𝒪
z
,
ɛ at the point z
∈
Ω. It leads to an asymptotic expansion of j in the following form:
(1)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
=
j
(
Ω
)
+
f
(
ɛ
)
g
(
z
)
+
o
(
f
(
ɛ
)
)
,
where f
(
ɛ
) is a scalar positive function going to zero with ɛ. This expression is called the topological asymptotic expansion and g is called the topological gradient.In order to minimizej, the best location to insert the obstacle 𝒪
z
,
ɛ in the fluid flow domain Ω is where g is the most negative. In fact if g
(
z
)
<
0, we have j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
<
j
(
Ω
).Starting with this observation, topological optimization algorithm can be constructed. The optimal domain is obtained using an iterative process building a sequence of geometries{
Ω
k
}
k
≥
0 with Ω
0
=
Ω. At the kth iteration, the topological gradient g
k is computed in Ω
k and the new geometry Ω
k
+
1 is obtained by inserting an obstacle 𝒪
k in the domain Ω
k; Ω
k
+
1
=
Ω
k
∖
𝒪
k
¯. The obstacle 𝒪
k is defined by a level set curve of g
k
(2)
𝒪
k
=
{
x
∈
Ω
k
,
such
that
g
k
(
x
)
≤
c
k
}
,
where c
k is chosen in such a way that the cost function j decreases as most as possible.To our knowledge, the topological sensitivity analysis for the nonstationary Navier-Stokes equations has not been studied so far. The most contributions have been focused on the stationary regime.The aim of this work is to extend the topological gradient method for the nonlinear unsteady flow environment. The main difficulty comes from the nonlinearity of the operator and the treatment of the associated adjoint problem. To overcome such a difficulty, we have used the discrete adjoint approach. This alternative takes a discretization of the Navier-Stokes equations, linearizes the discrete equations, and then uses the transpose of the linear operator to form the adjoint problem. The discrete adjoint approach has been developed by Elliott and Peraire [18], Mohammadi and Pironneau [19], and Nielsen and Anderson [20].The rest of this paper is organized as follows. In Section2, we present the topological optimization problem. In Section 3, we derive the topological asymptotic expansion for the discretized Navier-Stokes equations. The proposed numerical algorithm and some numerical results are described in Section 4. The paper ends by some concluding remarks.


## 2. Problem Statement

We consider a viscous and incompressible fluid flow in a bounded domainΩ of ℝ
3. The fluid flow is described by the nonstationary Navier-Stokes equations [21]. For given boundaries data (
u
d
,
g
) and a source term G (gravitational force), the velocity u
Ω and the pressure p
Ω satisfy the system
(3)
∂
u
Ω
∂
t
+
(
u
Ω
·
∇
)
u
Ω
-
ν
Δ
u
Ω
+
∇
p
Ω
=
G
in
Ω
×
(
0
,
T
)
,
∇
·
u
Ω
=
0
in
Ω
×
(
0
,
T
)
,
u
Ω
=
u
d
on
⁡
Γ
d
×
(
0
,
T
)
,
σ
(
u
Ω
,
p
Ω
)
·
n
=
g
on
⁡
Γ
n
×
(
0
,
T
)
,
u
Ω
(
·
,
0
)
=
u
0
in
Ω
,
where ν is the kinematic viscosity of the fluid, T is the computational time, and u
0 is an initial condition verifying ∇
·
u
0
=
0. Here Γ
d and Γ
n are two parts of the boundary ∂
Ω such that Γ
d
¯
∪
Γ
n
¯
=
∂
Ω and Γ
d
∩
Γ
n
=
∅. The stress distribution on Γ
n is defined by σ
(
u
,
p
)
·
n
=
g with σ
(
u
,
p
)
=
ν
(
∇
u
+
∇
u
T
)
-
p
I and n is the unit outward normal vector.


### 2.1. Shape Optimization Problem

We assume that the fluid flow domainΩ is a cavity having one inlet Γ
in and some outlets Γ
out
i (see Figure 1). The aim is to determine the optimal geometry Ω
p of the pipe bend domain minimizing the dissipated energy by the fluid flow. The considered problem can be formulated as follows:
(4)
(
O
)
{


Find
Ω
p
solution
to
min
⁡
D
∈
𝒟
a
d
j
(
D
)
,
where j is the cost function defined by
(5)
j
(
D
)
=
∫
0
T
‍
∫
Ω
‍
ν
|
∇
u
D
|
2
d
x
d
t
,
with u
D is the velocity field solution to the Navier-Stokes system (3) in the domain D.Figure 1
Cavity\Omega.Here𝒟
a
d is the set of admissible domains defined by
(6)
𝒟
a
d
=
{
D
⊂
Ω
suchthat
Γ
in
⊂
∂
Ω
∩
∂
D
,
Γ
out
i
⊂
∂
Ω
∩
∂
D
and
|
D
|
≤
V
desired
}
,
where |
D
| is the Lebesgue measure of D and V
desired denotes the target volume (weight).To solve the considered topology optimization problem(
O
), we propose an optimization approach based on the topological sensitivity analysis method. The optimal domain is constructed using an iterative process building a sequence of geometries {
D
k
}
k
≥
0 with D
0
=
Ω. At the kth iteration, the new geometry D
k
+
1 is obtained by inserting an obstacle 𝒪
k in the domain D
k; D
k
+
1
=
D
k
∖
𝒪
k
¯. The obstacle 𝒪
k is defined by a level set curve of a scalar function g
k
(7)
𝒪
k
=
{
x
∈
D
k
,
such
that
g
k
(
x
)
≤
c
k
}
,
where c
k is chosen in such a way that the cost function j decreases as most as possible. The function g
k is called the topological gradient. It is computed using a topological sensitivity analysis for the function j with respect to the creation of small geometry perturbation of the domain.


### 2.2. Optimal Location of a Small Obstacle

Let𝒪
z
,
ɛ be a small obstacle inserted in the fluid flow Ω. We assume that the obstacle has the form 𝒪
z
,
ɛ
=
z
+
ɛ
𝒪, where z
∈
Ω, ɛ
>
0, and 𝒪 is a given, fixed, and bounded domain of ℝ
3, containing the origin, whose boundary ∂
𝒪 is of 𝒞
1. In the presence of the obstacle, the function j is defined by
(8)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
=
∫
0
T
‍
∫
Ω
∖
𝒪
z
,
ɛ
¯
‍
ν
|
∇
u
ɛ
|
2
d
x
d
t
,
where u
ɛ is the velocity field solution to the Navier-Stokes system (9) in the perturbed domain Ω
z
,
ɛ
=
Ω
∖
𝒪
z
,
ɛ
¯:
(9)
∂
u
ɛ
∂
t
+
(
u
ɛ
·
∇
)
u
ɛ
-
ν
Δ
u
ɛ
+
∇
p
ɛ
=
𝒢
in
Ω
z
,
ɛ
×
(
0
,
T
)
,
∇
·
u
ɛ
=
0
in
Ω
z
,
ɛ
×
(
0
,
T
)
,
u
ɛ
=
u
d
on
⁡
Γ
d
×
(
0
,
T
)
,
σ
(
u
ɛ
,
p
ɛ
)
·
n
=
g
on
⁡
Γ
n
×
(
0
,
T
)
,
u
ɛ
=
0
on
⁡
∂
𝒪
z
,
ɛ
×
(
0
,
T
)
,
u
ɛ
(
·
,
0
)
=
u
0
in
Ω
z
,
ɛ
.Our aim is to determine the optimal location of the obstacle𝒪
z
,
ɛ in the fluid flow domain Ω minimizing the function j:
(10)
min
⁡
𝒪
z
,
ɛ
⊂
Ω
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
.
To solve this minimization problem, the topological sensitivity analysis method consists in studying the variation of the function j with respect to ɛ. It leads to an asymptotic expansion in the following form:
(11)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
=
j
(
Ω
)
+
f
(
ɛ
)
g
(
z
)
+
o
(
f
(
ɛ
)
)
,
where f
(
ɛ
) is a positive scalar function going to zero with ɛ. This expression is called the topological asymptotic expansion and g is called the topological gradient.From the asymptotic (11), one can observe that the best location of the obstacle 𝒪
z
,
ɛ in Ω minimizing j (solution to (10)) is given by
(12)
𝒪
z
⋆
,
ɛ
=
z
⋆
+
ɛ
𝒪
,
where z
⋆ is the location where g is the most negative. In fact if g
(
z
⋆
)
<
0 and g
(
z
⋆
)
≤
g
(
z
), for all z
∈
Ω, we have j
(
Ω
∖
𝒪
z
⋆
,
ɛ
¯
)
<
j
(
Ω
) and j
(
Ω
∖
𝒪
z
⋆
,
ɛ
¯
)
≤
j
(
Ω
∖
𝒪
z
,
ɛ
¯
), for all z
∈
Ω.Starting with this observation, the topological gradientg will be computed during the optimization process. It will be used to determine the location of the obstacle to be inserted at each iteration.To this end, we will derive a topological sensitivity analysis for the Navier-Stokes equations in the next section. The obtained results are valid for a large class of cost functionsj.


## 2.1. Shape Optimization Problem

We assume that the fluid flow domainΩ is a cavity having one inlet Γ
in and some outlets Γ
out
i (see Figure 1). The aim is to determine the optimal geometry Ω
p of the pipe bend domain minimizing the dissipated energy by the fluid flow. The considered problem can be formulated as follows:
(4)
(
O
)
{


Find
Ω
p
solution
to
min
⁡
D
∈
𝒟
a
d
j
(
D
)
,
where j is the cost function defined by
(5)
j
(
D
)
=
∫
0
T
‍
∫
Ω
‍
ν
|
∇
u
D
|
2
d
x
d
t
,
with u
D is the velocity field solution to the Navier-Stokes system (3) in the domain D.Figure 1
Cavity\Omega.Here𝒟
a
d is the set of admissible domains defined by
(6)
𝒟
a
d
=
{
D
⊂
Ω
suchthat
Γ
in
⊂
∂
Ω
∩
∂
D
,
Γ
out
i
⊂
∂
Ω
∩
∂
D
and
|
D
|
≤
V
desired
}
,
where |
D
| is the Lebesgue measure of D and V
desired denotes the target volume (weight).To solve the considered topology optimization problem(
O
), we propose an optimization approach based on the topological sensitivity analysis method. The optimal domain is constructed using an iterative process building a sequence of geometries {
D
k
}
k
≥
0 with D
0
=
Ω. At the kth iteration, the new geometry D
k
+
1 is obtained by inserting an obstacle 𝒪
k in the domain D
k; D
k
+
1
=
D
k
∖
𝒪
k
¯. The obstacle 𝒪
k is defined by a level set curve of a scalar function g
k
(7)
𝒪
k
=
{
x
∈
D
k
,
such
that
g
k
(
x
)
≤
c
k
}
,
where c
k is chosen in such a way that the cost function j decreases as most as possible. The function g
k is called the topological gradient. It is computed using a topological sensitivity analysis for the function j with respect to the creation of small geometry perturbation of the domain.


## 2.2. Optimal Location of a Small Obstacle

Let𝒪
z
,
ɛ be a small obstacle inserted in the fluid flow Ω. We assume that the obstacle has the form 𝒪
z
,
ɛ
=
z
+
ɛ
𝒪, where z
∈
Ω, ɛ
>
0, and 𝒪 is a given, fixed, and bounded domain of ℝ
3, containing the origin, whose boundary ∂
𝒪 is of 𝒞
1. In the presence of the obstacle, the function j is defined by
(8)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
=
∫
0
T
‍
∫
Ω
∖
𝒪
z
,
ɛ
¯
‍
ν
|
∇
u
ɛ
|
2
d
x
d
t
,
where u
ɛ is the velocity field solution to the Navier-Stokes system (9) in the perturbed domain Ω
z
,
ɛ
=
Ω
∖
𝒪
z
,
ɛ
¯:
(9)
∂
u
ɛ
∂
t
+
(
u
ɛ
·
∇
)
u
ɛ
-
ν
Δ
u
ɛ
+
∇
p
ɛ
=
𝒢
in
Ω
z
,
ɛ
×
(
0
,
T
)
,
∇
·
u
ɛ
=
0
in
Ω
z
,
ɛ
×
(
0
,
T
)
,
u
ɛ
=
u
d
on
⁡
Γ
d
×
(
0
,
T
)
,
σ
(
u
ɛ
,
p
ɛ
)
·
n
=
g
on
⁡
Γ
n
×
(
0
,
T
)
,
u
ɛ
=
0
on
⁡
∂
𝒪
z
,
ɛ
×
(
0
,
T
)
,
u
ɛ
(
·
,
0
)
=
u
0
in
Ω
z
,
ɛ
.Our aim is to determine the optimal location of the obstacle𝒪
z
,
ɛ in the fluid flow domain Ω minimizing the function j:
(10)
min
⁡
𝒪
z
,
ɛ
⊂
Ω
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
.
To solve this minimization problem, the topological sensitivity analysis method consists in studying the variation of the function j with respect to ɛ. It leads to an asymptotic expansion in the following form:
(11)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
=
j
(
Ω
)
+
f
(
ɛ
)
g
(
z
)
+
o
(
f
(
ɛ
)
)
,
where f
(
ɛ
) is a positive scalar function going to zero with ɛ. This expression is called the topological asymptotic expansion and g is called the topological gradient.From the asymptotic (11), one can observe that the best location of the obstacle 𝒪
z
,
ɛ in Ω minimizing j (solution to (10)) is given by
(12)
𝒪
z
⋆
,
ɛ
=
z
⋆
+
ɛ
𝒪
,
where z
⋆ is the location where g is the most negative. In fact if g
(
z
⋆
)
<
0 and g
(
z
⋆
)
≤
g
(
z
), for all z
∈
Ω, we have j
(
Ω
∖
𝒪
z
⋆
,
ɛ
¯
)
<
j
(
Ω
) and j
(
Ω
∖
𝒪
z
⋆
,
ɛ
¯
)
≤
j
(
Ω
∖
𝒪
z
,
ɛ
¯
), for all z
∈
Ω.Starting with this observation, the topological gradientg will be computed during the optimization process. It will be used to determine the location of the obstacle to be inserted at each iteration.To this end, we will derive a topological sensitivity analysis for the Navier-Stokes equations in the next section. The obtained results are valid for a large class of cost functionsj.


## 3. Topological Sensitivity Analysis

We start our analysis by the time discretization of the Navier-Stokes problem. It leads to solving steady state generalized Stokes equations at each time step. The topological sensitivity analysis for the Navier-Stokes equations is derived in Section3.3.


### 3.1. The Time Discretization

We remark that the convective term in the first equation of system (9) coincides with the total derivative,
(13)
d
u
ɛ
d
t
=
∂
u
ɛ
∂
t
+
(
u
ɛ
·
∇
)
u
ɛ
.
Then, if Δ
t is the time step and u
ɛ
n is the approximated velocity at time t
n
=
n
Δ
t, using the characteristics method [22], we have the following approximation:
(14)
d
u
ɛ
d
t
(
x
,
t
n
+
1
)
=
u
ɛ
n
+
1
(
x
)
-
u
ɛ
n
∘
χ
n
(
x
)
Δ
t
,
where χ
n
(
x
)
=
X
n
(
t
n
+
1
;
x
) describes the position at time t
n
+
1
=
(
n
+
1
)
Δ
t of the particle of fluid which is at point x at time t
n. It is the solution of
(15)
d
X
n
d
τ
(
τ
;
x
)
=
u
ɛ
n
(
X
n
(
τ
;
x
)
,
τ
)
,
t
n
<
τ
<
t
n
+
1
X
(
t
n
;
x
)
=
x
.
Using (15), the time discretization of (9) reads
(16)
c
u
ɛ
n
+
1
-
ν
Δ
u
ɛ
n
+
1
+
∇
p
ɛ
n
+
1
=
F
n
+
1
in
Ω
z
,
ɛ
,
∇
·
u
ɛ
n
+
1
=
0
in
Ω
z
,
ɛ
,
u
ɛ
n
+
1
=
u
d
on
⁡
Γ
d
,
σ
(
u
ɛ
n
+
1
,
p
ɛ
n
+
1
)
·
n
=
g
on
⁡
Γ
n
,
u
ɛ
n
+
1
=
0
on
⁡
∂
𝒪
z
,
ɛ
,
where c
=
1
/
Δ
t, F
n
+
1
=
(
1
/
Δ
t
)
u
ɛ
n
o
χ
n
+
𝒢
n
+
1, u
ɛ
n
+
1, and p
ɛ
n
+
1 are the approximations of u
ɛ and p
ɛ on time t
n
+
1.Then, at each time step, we have to solve a steady state generalized Stokes problem having the following generic form:(17)
c
u
ɛ
-
ν
Δ
u
ɛ
+
∇
p
ɛ
=
F
in
Ω
z
,
ɛ
,
∇
·
u
ɛ
=
0
in
Ω
z
,
ɛ
,
u
ɛ
=
u
d
on
⁡
Γ
d
,
σ
(
u
ɛ
,
p
ɛ
)
·
n
=
g
on
⁡
Γ
n
,
u
ɛ
=
0
on
⁡
∂
𝒪
z
,
ɛ
.
For u
ɛ
∈
{
v
∈
H
1
(
Ω
z
,
ɛ
)
,
∇
·
v
=
0
in
Ω
z
,
ɛ
)
}, we show in [10] that problem (17) has a unique solution.


### 3.2. Topological Sensitivity Analysis for the Generalized Stokes Equations

In this section, we give the topological sensitivity analysis for the generalized Stokes equations when creating a small hole𝒪
z
,
ɛ inside the domain with a homogeneous Dirichlet condition on the boundary ∂
𝒪
ɛ. We recall here the main results of this case. For more details, one can consult [10] or [15]. The presented results are valid for all functional J
ɛ
:
H
1
(
Ω
z
,
ɛ
)
3
→
ℝ satisfying the following assumptions.Hypothesis 1.
(i)
The functionJ
0 is differentiable with respect to u, and we denote by D
J
0
(
u
) its derivative.
(ii)
There exists a real numberδ
J such that for all ɛ
≥
0
(18)
J
ɛ
(
u
ɛ
)
-
J
0
(
u
0
)
=
D
J
0
(
u
0
)
(
u
ɛ
-
u
0
)
+
ɛ
δ
J
+
o
(
ɛ
)
.Theorem 1 (see [10, 15]).
Under the assumptions of Hypothesis1, the function J
ɛ has the asymptotic expansion
(19)
J
ɛ
(
u
ɛ
)
=
J
0
(
u
0
)
+
ɛ
[
(
-
∫
∂
𝒪
‍
η
(
y
)
d
s
(
y
)
)
·
v
0
(
z
)
+
δ
J
(
z
)
]
+
o
(
ɛ
)
.
The function v
0 is the solution to the adjoint problem
(20)
c
v
0
-
ν
Δ
v
0
+
∇
q
0
=
-
D
J
0
(
u
0
)
i
n
Ω
,
∇
·
v
0
=
0
i
n
Ω
,
v
0
=
0
o
n
Γ
d
,
σ
(
v
0
,
q
0
)
·
n
=
0
o
n
Γ
n
.
The function η
∈
H
-
1
/
2
(
∂
𝒪
)
3 is solution to the following boundary integral equation (for more details, one can see [15, 23]):
(21)
∫
∂
𝒪
‍
E
(
x
-
y
)
η
(
y
)
d
s
(
y
)
=
-
u
0
(
z
)
,
∀
x
∈
∂
𝒪
,
where (E, P) is the fundamental solution to the Stokes equations in ℝ
3
(22)
E
(
y
)
=
(
1
8
π
ν
r
)
(
I
+
e
r
e
r
T
)
,
P
(
y
)
=
y
4
π
r
3
,
where r
=
∥
y
∥, e
r
=
y
/
r, and e
r
T is the transposed vector of e
r.
In the particular case where𝒪 is the unit ball B
(
0,1
), we have
(23)
∫
∂
ω
‍
E
(
x
-
y
)
d
s
(
y
)
=
2
3
ν
I
,
∀
x
∈
∂
𝒪
.
Hence, the density η is given explicitly η
(
y
)
=
-
(
3
ν
/
2
)
u
0
(
z
), ∀
y
∈
∂
𝒪.Corollary 2.
Letz
∈
Ω and 𝒪
=
B
(
0,1
). Under the hypotheses of Theorem 1, one has
(24)
J
ɛ
(
u
ɛ
)
=
J
0
(
u
0
)
+
ɛ
[
6
π
ν
u
0
(
z
)
·
v
0
(
z
)
+
δ
J
(
z
)
]
+
o
(
ɛ
)
.
The expression of the term δ
J depends on the considered cost function J
ɛ. In the following proposition, we consider the seminorm and we calculate its variation δ
J.Proposition 3 (see [10]).
Letw
d
∈
H
2
(
Ω
) be a given wanted velocity field.
The cost functionJ
ɛ
(
u
)
=
ν
∫
Ω
∖
𝒪
z
,
ɛ
¯
‍
|
∇
u
-
∇
w
d
|
2
d
x satisfies the assumptions of Hypothesis 1 with
(25)
D
J
0
(
w
)
=
2
∫
Ω
‍
∇
(
u
0
-
w
d
)
∇
w
d
x
a
n
d
δ
J
(
z
)
=
(
-
∫
∂
𝒪
‍
η
(
y
)
d
s
(
y
)
)
·
v
0
(
z
)
.
If 𝒪 is the unit ball B
(
0,1
), we have δ
J
(
z
)
=
6
π
ν
|
u
0
(
z
)
|
2, for all z
∈
Ω.


### 3.3. Topological Sensitivity Analysis for the Discretized Navier-Stokes Equations

Consider a shape functionj of the form
(26)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
=
𝒥
ɛ
(
u
ɛ
)
,
with 𝒥
ɛ
(
u
ɛ
)
=
∫
0
T
‍
J
ɛ
(
u
ɛ
)
d
t for all ɛ
≥
0, where u
ɛ
∈
L
2
(
0
,
T
;
H
1
(
Ω
z
,
ɛ
)
3
) is the solution to (9) and the functional J
ɛ
:
H
1
(
Ω
z
,
ɛ
)
3
→
ℝ satisfies the following assumptions.Hypothesis 2.
(i)
J
0 is differentiable with respect to u, with D
J
0
(
u
0
(
·
,
t
)
)
∈
H
-
1
(
Ω
)
3 for almost all t
∈
(
0
,
T
).
(ii)
For allɛ
≥
0, J
ɛ
(
u
)
∈
L
1
(
0
,
T
), for all u
∈
H
1
(
Ω
z
,
ɛ
)
3.
(iii)
There exist a real numberδ
J, independent of ɛ, such that for all ɛ
≥
0
(27)
𝒥
ɛ
(
u
ɛ
)
-
𝒥
0
(
u
0
)
=
∫
0
T
‍
〈
D
J
0
(
u
0
)
,
u
ɛ
-
u
0
〉
H
-
1
(
Ω
)
×
H
0
1
(
Ω
)
d
t
+
ɛ
δ
J
+
o
(
ɛ
)
.
In this section we consider the nonstationary Navier-Stokes equations and we compute the variation of the cost functionj with respect to the insertion of a small obstacle in the fluid flow domain. The basic idea is to use the discretized formulation (16) and apply the established results for the generalized Stokes system. We have
(28)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
-
j
(
Ω
)
=
∫
0
T
‍
J
ɛ
(
u
ɛ
)
d
t
-
∫
0
T
‍
J
0
(
u
0
)
d
t
=
∑
n
=
1
N
T
‍
∫
t
n
-
1
t
n
‍
(
J
ɛ
(
u
ɛ
)
-
J
0
(
u
0
)
)
d
t
,
where N
T
∈
ℕ
* and t
n
=
n
Δ
t, with Δ
t
=
T
/
N
T being the time step. From the discrete formulation of the Navier-Stokes equations (16), we obtain the following approximation:
(29)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
-
j
(
Ω
)
=
Δ
t
∑
n
=
1
N
T
‍
[
J
ɛ
(
u
ɛ
n
)
-
J
0
(
u
0
n
)
]
.
Using the sensitivity analysis for the generalized Stokes equations, presented in the last section, we deduce
(30)
J
ɛ
(
u
ɛ
n
)
-
J
0
(
u
0
n
)
=
f
(
ɛ
)
Δ
t
δ
j
n
(
z
)
+
o
(
f
(
ɛ
)
)
,
with
(31)
δ
j
n
(
z
)
=
(
-
∫
∂
𝒪
‍
η
n
(
y
)
d
s
(
y
)
)
·
v
0
n
(
z
)
+
δ
J
n
(
z
)
∀
z
∈
Ω
,
where u
0
n is the solution to the discretized Navier-Stokes equations at time t
n
=
n
Δ
t
(32)
c
u
0
n
-
ν
Δ
u
0
n
+
∇
p
0
n
=
F
n
in
Ω
,
∇
·
u
0
n
=
0
in
Ω
,
u
0
n
=
u
d
on
⁡
Γ
d
,
σ
(
u
0
,
p
0
)
·
n
=
g
on
⁡
Γ
n
.
The field v
0
n is the solution to the associated adjoint problem. It is computed in reverse time from the final flow solution
(33)
c
v
0
n
-
ν
Δ
v
0
n
+
∇
q
0
n
=
-
D
J
0
(
u
0
N
T
-
n
)
in
Ω
,
∇
·
v
0
n
=
0
in
Ω
,
v
0
n
=
0
on
⁡
Γ
d
,
σ
(
v
0
,
q
0
)
·
n
=
0
on
⁡
Γ
n
,
where u
0
N
T
-
n is the approximation of the velocity field u
0
(
·
,
T
-
t
) at time T
-
t
n.
The functionη
n
∈
H
-
1
/
2
(
∂
𝒪
)
d is solution to the boundary integral equation
(34)
∫
∂
𝒪
‍
E
(
x
-
y
)
η
n
(
y
)
d
s
(
y
)
=
-
u
0
n
(
z
)
,
∀
x
∈
∂
𝒪
.
Then, we derive the following results for the Navier-Stokes equation.Theorem 4.
IfJ
ɛ satisfies the assumptions of Hypothesis 2, the function j admits the following asymptotic expansion:
(35)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
-
j
(
Ω
)
=
f
(
ɛ
)
Δ
t
g
(
z
)
+
o
(
f
(
ɛ
)
)
,
where g is the topological gradient given by
(36)
g
(
z
)
=
∑
n
=
1
N
T
‍
δ
j
n
(
z
)
,
∀
z
∈
Ω
.
If 𝒪 is the unit ball B
(
0,1
), we have δ
j
n
(
z
)
=
6
π
ν
u
0
n
(
z
)
·
v
0
n
(
z
)
+
δ
J
n
(
z
).
Then, we deduce the following corollary.Corollary 5.
Letz
∈
Ω and 𝒪
=
B
(
0,1
). Under the hypotheses of Theorem 4, one has
(37)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
-
j
(
Ω
)
=
ɛ
Δ
t
∑
n
=
1
N
T
‍
[
6
π
ν
u
0
n
(
z
)
·
v
0
n
(
z
)
+
δ
J
n
(
z
)
]
+
o
(
ɛ
)
.


## 3.1. The Time Discretization

We remark that the convective term in the first equation of system (9) coincides with the total derivative,
(13)
d
u
ɛ
d
t
=
∂
u
ɛ
∂
t
+
(
u
ɛ
·
∇
)
u
ɛ
.
Then, if Δ
t is the time step and u
ɛ
n is the approximated velocity at time t
n
=
n
Δ
t, using the characteristics method [22], we have the following approximation:
(14)
d
u
ɛ
d
t
(
x
,
t
n
+
1
)
=
u
ɛ
n
+
1
(
x
)
-
u
ɛ
n
∘
χ
n
(
x
)
Δ
t
,
where χ
n
(
x
)
=
X
n
(
t
n
+
1
;
x
) describes the position at time t
n
+
1
=
(
n
+
1
)
Δ
t of the particle of fluid which is at point x at time t
n. It is the solution of
(15)
d
X
n
d
τ
(
τ
;
x
)
=
u
ɛ
n
(
X
n
(
τ
;
x
)
,
τ
)
,
t
n
<
τ
<
t
n
+
1
X
(
t
n
;
x
)
=
x
.
Using (15), the time discretization of (9) reads
(16)
c
u
ɛ
n
+
1
-
ν
Δ
u
ɛ
n
+
1
+
∇
p
ɛ
n
+
1
=
F
n
+
1
in
Ω
z
,
ɛ
,
∇
·
u
ɛ
n
+
1
=
0
in
Ω
z
,
ɛ
,
u
ɛ
n
+
1
=
u
d
on
⁡
Γ
d
,
σ
(
u
ɛ
n
+
1
,
p
ɛ
n
+
1
)
·
n
=
g
on
⁡
Γ
n
,
u
ɛ
n
+
1
=
0
on
⁡
∂
𝒪
z
,
ɛ
,
where c
=
1
/
Δ
t, F
n
+
1
=
(
1
/
Δ
t
)
u
ɛ
n
o
χ
n
+
𝒢
n
+
1, u
ɛ
n
+
1, and p
ɛ
n
+
1 are the approximations of u
ɛ and p
ɛ on time t
n
+
1.Then, at each time step, we have to solve a steady state generalized Stokes problem having the following generic form:(17)
c
u
ɛ
-
ν
Δ
u
ɛ
+
∇
p
ɛ
=
F
in
Ω
z
,
ɛ
,
∇
·
u
ɛ
=
0
in
Ω
z
,
ɛ
,
u
ɛ
=
u
d
on
⁡
Γ
d
,
σ
(
u
ɛ
,
p
ɛ
)
·
n
=
g
on
⁡
Γ
n
,
u
ɛ
=
0
on
⁡
∂
𝒪
z
,
ɛ
.
For u
ɛ
∈
{
v
∈
H
1
(
Ω
z
,
ɛ
)
,
∇
·
v
=
0
in
Ω
z
,
ɛ
)
}, we show in [10] that problem (17) has a unique solution.


## 3.2. Topological Sensitivity Analysis for the Generalized Stokes Equations

In this section, we give the topological sensitivity analysis for the generalized Stokes equations when creating a small hole𝒪
z
,
ɛ inside the domain with a homogeneous Dirichlet condition on the boundary ∂
𝒪
ɛ. We recall here the main results of this case. For more details, one can consult [10] or [15]. The presented results are valid for all functional J
ɛ
:
H
1
(
Ω
z
,
ɛ
)
3
→
ℝ satisfying the following assumptions.Hypothesis 1.
(i)
The functionJ
0 is differentiable with respect to u, and we denote by D
J
0
(
u
) its derivative.
(ii)
There exists a real numberδ
J such that for all ɛ
≥
0
(18)
J
ɛ
(
u
ɛ
)
-
J
0
(
u
0
)
=
D
J
0
(
u
0
)
(
u
ɛ
-
u
0
)
+
ɛ
δ
J
+
o
(
ɛ
)
.Theorem 1 (see [10, 15]).
Under the assumptions of Hypothesis1, the function J
ɛ has the asymptotic expansion
(19)
J
ɛ
(
u
ɛ
)
=
J
0
(
u
0
)
+
ɛ
[
(
-
∫
∂
𝒪
‍
η
(
y
)
d
s
(
y
)
)
·
v
0
(
z
)
+
δ
J
(
z
)
]
+
o
(
ɛ
)
.
The function v
0 is the solution to the adjoint problem
(20)
c
v
0
-
ν
Δ
v
0
+
∇
q
0
=
-
D
J
0
(
u
0
)
i
n
Ω
,
∇
·
v
0
=
0
i
n
Ω
,
v
0
=
0
o
n
Γ
d
,
σ
(
v
0
,
q
0
)
·
n
=
0
o
n
Γ
n
.
The function η
∈
H
-
1
/
2
(
∂
𝒪
)
3 is solution to the following boundary integral equation (for more details, one can see [15, 23]):
(21)
∫
∂
𝒪
‍
E
(
x
-
y
)
η
(
y
)
d
s
(
y
)
=
-
u
0
(
z
)
,
∀
x
∈
∂
𝒪
,
where (E, P) is the fundamental solution to the Stokes equations in ℝ
3
(22)
E
(
y
)
=
(
1
8
π
ν
r
)
(
I
+
e
r
e
r
T
)
,
P
(
y
)
=
y
4
π
r
3
,
where r
=
∥
y
∥, e
r
=
y
/
r, and e
r
T is the transposed vector of e
r.
In the particular case where𝒪 is the unit ball B
(
0,1
), we have
(23)
∫
∂
ω
‍
E
(
x
-
y
)
d
s
(
y
)
=
2
3
ν
I
,
∀
x
∈
∂
𝒪
.
Hence, the density η is given explicitly η
(
y
)
=
-
(
3
ν
/
2
)
u
0
(
z
), ∀
y
∈
∂
𝒪.Corollary 2.
Letz
∈
Ω and 𝒪
=
B
(
0,1
). Under the hypotheses of Theorem 1, one has
(24)
J
ɛ
(
u
ɛ
)
=
J
0
(
u
0
)
+
ɛ
[
6
π
ν
u
0
(
z
)
·
v
0
(
z
)
+
δ
J
(
z
)
]
+
o
(
ɛ
)
.
The expression of the term δ
J depends on the considered cost function J
ɛ. In the following proposition, we consider the seminorm and we calculate its variation δ
J.Proposition 3 (see [10]).
Letw
d
∈
H
2
(
Ω
) be a given wanted velocity field.
The cost functionJ
ɛ
(
u
)
=
ν
∫
Ω
∖
𝒪
z
,
ɛ
¯
‍
|
∇
u
-
∇
w
d
|
2
d
x satisfies the assumptions of Hypothesis 1 with
(25)
D
J
0
(
w
)
=
2
∫
Ω
‍
∇
(
u
0
-
w
d
)
∇
w
d
x
a
n
d
δ
J
(
z
)
=
(
-
∫
∂
𝒪
‍
η
(
y
)
d
s
(
y
)
)
·
v
0
(
z
)
.
If 𝒪 is the unit ball B
(
0,1
), we have δ
J
(
z
)
=
6
π
ν
|
u
0
(
z
)
|
2, for all z
∈
Ω.


## 3.3. Topological Sensitivity Analysis for the Discretized Navier-Stokes Equations

Consider a shape functionj of the form
(26)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
=
𝒥
ɛ
(
u
ɛ
)
,
with 𝒥
ɛ
(
u
ɛ
)
=
∫
0
T
‍
J
ɛ
(
u
ɛ
)
d
t for all ɛ
≥
0, where u
ɛ
∈
L
2
(
0
,
T
;
H
1
(
Ω
z
,
ɛ
)
3
) is the solution to (9) and the functional J
ɛ
:
H
1
(
Ω
z
,
ɛ
)
3
→
ℝ satisfies the following assumptions.Hypothesis 2.
(i)
J
0 is differentiable with respect to u, with D
J
0
(
u
0
(
·
,
t
)
)
∈
H
-
1
(
Ω
)
3 for almost all t
∈
(
0
,
T
).
(ii)
For allɛ
≥
0, J
ɛ
(
u
)
∈
L
1
(
0
,
T
), for all u
∈
H
1
(
Ω
z
,
ɛ
)
3.
(iii)
There exist a real numberδ
J, independent of ɛ, such that for all ɛ
≥
0
(27)
𝒥
ɛ
(
u
ɛ
)
-
𝒥
0
(
u
0
)
=
∫
0
T
‍
〈
D
J
0
(
u
0
)
,
u
ɛ
-
u
0
〉
H
-
1
(
Ω
)
×
H
0
1
(
Ω
)
d
t
+
ɛ
δ
J
+
o
(
ɛ
)
.
In this section we consider the nonstationary Navier-Stokes equations and we compute the variation of the cost functionj with respect to the insertion of a small obstacle in the fluid flow domain. The basic idea is to use the discretized formulation (16) and apply the established results for the generalized Stokes system. We have
(28)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
-
j
(
Ω
)
=
∫
0
T
‍
J
ɛ
(
u
ɛ
)
d
t
-
∫
0
T
‍
J
0
(
u
0
)
d
t
=
∑
n
=
1
N
T
‍
∫
t
n
-
1
t
n
‍
(
J
ɛ
(
u
ɛ
)
-
J
0
(
u
0
)
)
d
t
,
where N
T
∈
ℕ
* and t
n
=
n
Δ
t, with Δ
t
=
T
/
N
T being the time step. From the discrete formulation of the Navier-Stokes equations (16), we obtain the following approximation:
(29)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
-
j
(
Ω
)
=
Δ
t
∑
n
=
1
N
T
‍
[
J
ɛ
(
u
ɛ
n
)
-
J
0
(
u
0
n
)
]
.
Using the sensitivity analysis for the generalized Stokes equations, presented in the last section, we deduce
(30)
J
ɛ
(
u
ɛ
n
)
-
J
0
(
u
0
n
)
=
f
(
ɛ
)
Δ
t
δ
j
n
(
z
)
+
o
(
f
(
ɛ
)
)
,
with
(31)
δ
j
n
(
z
)
=
(
-
∫
∂
𝒪
‍
η
n
(
y
)
d
s
(
y
)
)
·
v
0
n
(
z
)
+
δ
J
n
(
z
)
∀
z
∈
Ω
,
where u
0
n is the solution to the discretized Navier-Stokes equations at time t
n
=
n
Δ
t
(32)
c
u
0
n
-
ν
Δ
u
0
n
+
∇
p
0
n
=
F
n
in
Ω
,
∇
·
u
0
n
=
0
in
Ω
,
u
0
n
=
u
d
on
⁡
Γ
d
,
σ
(
u
0
,
p
0
)
·
n
=
g
on
⁡
Γ
n
.
The field v
0
n is the solution to the associated adjoint problem. It is computed in reverse time from the final flow solution
(33)
c
v
0
n
-
ν
Δ
v
0
n
+
∇
q
0
n
=
-
D
J
0
(
u
0
N
T
-
n
)
in
Ω
,
∇
·
v
0
n
=
0
in
Ω
,
v
0
n
=
0
on
⁡
Γ
d
,
σ
(
v
0
,
q
0
)
·
n
=
0
on
⁡
Γ
n
,
where u
0
N
T
-
n is the approximation of the velocity field u
0
(
·
,
T
-
t
) at time T
-
t
n.
The functionη
n
∈
H
-
1
/
2
(
∂
𝒪
)
d is solution to the boundary integral equation
(34)
∫
∂
𝒪
‍
E
(
x
-
y
)
η
n
(
y
)
d
s
(
y
)
=
-
u
0
n
(
z
)
,
∀
x
∈
∂
𝒪
.
Then, we derive the following results for the Navier-Stokes equation.Theorem 4.
IfJ
ɛ satisfies the assumptions of Hypothesis 2, the function j admits the following asymptotic expansion:
(35)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
-
j
(
Ω
)
=
f
(
ɛ
)
Δ
t
g
(
z
)
+
o
(
f
(
ɛ
)
)
,
where g is the topological gradient given by
(36)
g
(
z
)
=
∑
n
=
1
N
T
‍
δ
j
n
(
z
)
,
∀
z
∈
Ω
.
If 𝒪 is the unit ball B
(
0,1
), we have δ
j
n
(
z
)
=
6
π
ν
u
0
n
(
z
)
·
v
0
n
(
z
)
+
δ
J
n
(
z
).
Then, we deduce the following corollary.Corollary 5.
Letz
∈
Ω and 𝒪
=
B
(
0,1
). Under the hypotheses of Theorem 4, one has
(37)
j
(
Ω
∖
𝒪
z
,
ɛ
¯
)
-
j
(
Ω
)
=
ɛ
Δ
t
∑
n
=
1
N
T
‍
[
6
π
ν
u
0
n
(
z
)
·
v
0
n
(
z
)
+
δ
J
n
(
z
)
]
+
o
(
ɛ
)
.


## 4. Algorithm and Numerical Results

This section is devoted to some numerical investigations for our shape optimization problem(
O
) (see Section 2). The optimization algorithm is based on the topological sensitivity analysis obtained in Section 3. From Corollary 5 and Proposition 3, we deduce that the function j
(
D
)
=
∫
0
T
‍
∫
Ω
‍
|
∇
u
|
2
d
x
d
t admits the following asymptotic expansion:
(38)
j
(
D
∖
𝒪
z
,
ɛ
¯
)
-
j
(
D
)
=
ɛ
Δ
t
g
(
z
)
+
o
(
ɛ
)
,
where the topological gradient g is given by
(39)
g
(
z
)
=
6
π
ν
∑
n
=
1
N
T
‍
(
u
0
n
(
z
)
·
v
0
n
(
z
)
+
|
u
0
n
(
z
)
|
2
)
,
∀
z
∈
Ω
.
The expression of g is deduced from Corollary 5 and Proposition 3.


### 4.1. The Topological Optimization Algorithm

As already mentioned, the optimal domain is obtained using an iterative process building a sequence of geometries(
Ω
k
)
k
≥
0 with Ω
0
=
Ω. At the kth, iteration the topological gradient is denoted by g
k and the new geometry Ω
k
+
1 is obtained by inserting an obstacle 𝒪
k in the domain Ω
k; Ω
k
+
1
=
Ω
k
∖
𝒪
k
¯. The obstacle 𝒪
k is defined by a level set curve of the topological gradient g
k.The Algorithm. Topology optimization with volume constraint.(i)
Initialization: chooseΩ
0
=
Ω, and set k
=
0.
(ii)
Repeat until|
Ω
k
|
≤
V
desired:
(a)
computeu
k
=
(
u
k
n
)
1
≤
n
≤
N
T the solution to the Navier-Stokes equations (32) in Ω
k,
(b)
computev
k
=
(
v
k
n
)
1
≤
n
≤
N
T the solution to the associated adjoint problem (33) in Ω
k,
(c)
compute the topological gradientg
k
(
z
), ∀
z
∈
Ω
k,
(d)
determine the obstacle𝒪
k
=
{
x
∈
Ω
k
,
such
that
g
k
(
x
)
≤
c
k
<
0
},
(e)
get the new domainΩ
k
+
1
=
Ω
k
∖
𝒪
k
¯
(f)
k
←
k
+
1.The functiong
k is computed piecewise constant over elements. The constant c
k determines the volume of the obstacle 𝒪
k to be inserted. In practice, c
k is chosen in such a way that the obstacle volume |
𝒪
k
| is less or equal to 10
% of the current domain volume |
Ω
k
|; that is, |
𝒪
k
|
/
|
Ω
k
|
≤
0.1. More precisely, 𝒪
k is inserted in the zone Ω
k
-
=
{
x
∈
Ω
k
,
g
k
(
x
)
<
0
} where g
k is negative. The parameter c
k is chosen such that |
𝒪
k
|
=
min
⁡
{
|
Ω
k
-
|
,
0.1
|
Ω
k
|
}. The condition 𝒪
k
⊂
Ω
k
- ensured the decrease of the design function j
(
Ω
k
+
1
)
<
j
(
Ω
k
).


#### 4.1.1. A Cavity with One Inlet and Two Outlets

We consider a cavity with one inlet and two outlets having the same section (see Figure2(a)) where Γ
in is a disc of center z
in
=
(
0,0.5,0.5
) and radius r
in
=
0.1, Γ
out
1 is a disc of center z
out
1
=
(
0.75,0.5,0
) and radius r
out
=
0.05, and Γ
out
2 is a disc of center z
out
2
=
(
0.75,0.5,1
) and radius r
out
=
0.05. The results of this case are given in Figures 3 and 4. The optimal geometry is presented in Figure 3(a). It is obtained for V
desired
=
0.15
|
Ω
|. A vertical cut of the obtained velocity is shown in Figure 3(b). The variation of 𝒥 is illustrated in Figure 3(c). We present in Figure 4 some geometries obtained during the optimization process.The initial domain and Figureelocity field.
(a)
The initial domain
(b)
The initial velocity fieldFirst case: cavity with one inlet and two outlets.
(a)
The optimal design
(b)
The obtained velocity field
(c)
Variation of the function𝒥The obtained geometries during the optimization process.
(a)
Iteration 2
(b)
Iteration 3
(c)
Iteration 6
(d)
Iteration 8
(e)
Iteration 13
(f)
Iteration 19


#### 4.1.2. A Cavity with One Inlet and Three Outlets

In this case, we use a cavity with one inletΓ
in and three outlets Γ
out
1, Γ
out
2, and Γ
out
3 (see Figure 5(a)). The inlet and the outlets are defined as(i)
Γ
in is a disc of center (0,0.5,0.2) and diameter equal to 0.2;
(ii)
Γ
out
1 is a disc of center (1,0.5,0.4) and diameter equal to 0.066;
(iii)
Γ
out
2 is a disc of center (1,0.5,0.6) and diameter equal to 0.066;
(iv)
Γ
out
3 is a disc of center (1,0.5,0.8) and diameter equal to 0.066.The initial domain and velocity field.
(a)
The initial domain
(b)
The initial velocity fieldThe initial domain and the initial velocity field are given in Figure5. The results of this example are illustrated in Figures 6 and 7. The optimal domain is computed for V
desired
=
0.2
|
Ω
|. It is presented in Figure 6(a). A vertical cut (at y
=
0.5) of the velocity is shown in Figure 6(b). Figure 7 shows some geometries obtained during the optimization process.Second case: a cavity with one inlet and three outlets.
(a)
The optimal shape design
(b)
The obtained velocity fieldThe obtained geometries during the optimization process.
(a)
Iteration 2
(b)
Iteration 3
(c)
Iteration 5
(d)
Iteration 8
(e)
Iteration 11
(f)
Iteration 17


## 4.1. The Topological Optimization Algorithm

As already mentioned, the optimal domain is obtained using an iterative process building a sequence of geometries(
Ω
k
)
k
≥
0 with Ω
0
=
Ω. At the kth, iteration the topological gradient is denoted by g
k and the new geometry Ω
k
+
1 is obtained by inserting an obstacle 𝒪
k in the domain Ω
k; Ω
k
+
1
=
Ω
k
∖
𝒪
k
¯. The obstacle 𝒪
k is defined by a level set curve of the topological gradient g
k.The Algorithm. Topology optimization with volume constraint.(i)
Initialization: chooseΩ
0
=
Ω, and set k
=
0.
(ii)
Repeat until|
Ω
k
|
≤
V
desired:
(a)
computeu
k
=
(
u
k
n
)
1
≤
n
≤
N
T the solution to the Navier-Stokes equations (32) in Ω
k,
(b)
computev
k
=
(
v
k
n
)
1
≤
n
≤
N
T the solution to the associated adjoint problem (33) in Ω
k,
(c)
compute the topological gradientg
k
(
z
), ∀
z
∈
Ω
k,
(d)
determine the obstacle𝒪
k
=
{
x
∈
Ω
k
,
such
that
g
k
(
x
)
≤
c
k
<
0
},
(e)
get the new domainΩ
k
+
1
=
Ω
k
∖
𝒪
k
¯
(f)
k
←
k
+
1.The functiong
k is computed piecewise constant over elements. The constant c
k determines the volume of the obstacle 𝒪
k to be inserted. In practice, c
k is chosen in such a way that the obstacle volume |
𝒪
k
| is less or equal to 10
% of the current domain volume |
Ω
k
|; that is, |
𝒪
k
|
/
|
Ω
k
|
≤
0.1. More precisely, 𝒪
k is inserted in the zone Ω
k
-
=
{
x
∈
Ω
k
,
g
k
(
x
)
<
0
} where g
k is negative. The parameter c
k is chosen such that |
𝒪
k
|
=
min
⁡
{
|
Ω
k
-
|
,
0.1
|
Ω
k
|
}. The condition 𝒪
k
⊂
Ω
k
- ensured the decrease of the design function j
(
Ω
k
+
1
)
<
j
(
Ω
k
).


### 4.1.1. A Cavity with One Inlet and Two Outlets

We consider a cavity with one inlet and two outlets having the same section (see Figure2(a)) where Γ
in is a disc of center z
in
=
(
0,0.5,0.5
) and radius r
in
=
0.1, Γ
out
1 is a disc of center z
out
1
=
(
0.75,0.5,0
) and radius r
out
=
0.05, and Γ
out
2 is a disc of center z
out
2
=
(
0.75,0.5,1
) and radius r
out
=
0.05. The results of this case are given in Figures 3 and 4. The optimal geometry is presented in Figure 3(a). It is obtained for V
desired
=
0.15
|
Ω
|. A vertical cut of the obtained velocity is shown in Figure 3(b). The variation of 𝒥 is illustrated in Figure 3(c). We present in Figure 4 some geometries obtained during the optimization process.The initial domain and Figureelocity field.
(a)
The initial domain
(b)
The initial velocity fieldFirst case: cavity with one inlet and two outlets.
(a)
The optimal design
(b)
The obtained velocity field
(c)
Variation of the function𝒥The obtained geometries during the optimization process.
(a)
Iteration 2
(b)
Iteration 3
(c)
Iteration 6
(d)
Iteration 8
(e)
Iteration 13
(f)
Iteration 19


### 4.1.2. A Cavity with One Inlet and Three Outlets

In this case, we use a cavity with one inletΓ
in and three outlets Γ
out
1, Γ
out
2, and Γ
out
3 (see Figure 5(a)). The inlet and the outlets are defined as(i)
Γ
in is a disc of center (0,0.5,0.2) and diameter equal to 0.2;
(ii)
Γ
out
1 is a disc of center (1,0.5,0.4) and diameter equal to 0.066;
(iii)
Γ
out
2 is a disc of center (1,0.5,0.6) and diameter equal to 0.066;
(iv)
Γ
out
3 is a disc of center (1,0.5,0.8) and diameter equal to 0.066.The initial domain and velocity field.
(a)
The initial domain
(b)
The initial velocity fieldThe initial domain and the initial velocity field are given in Figure5. The results of this example are illustrated in Figures 6 and 7. The optimal domain is computed for V
desired
=
0.2
|
Ω
|. It is presented in Figure 6(a). A vertical cut (at y
=
0.5) of the velocity is shown in Figure 6(b). Figure 7 shows some geometries obtained during the optimization process.Second case: a cavity with one inlet and three outlets.
(a)
The optimal shape design
(b)
The obtained velocity fieldThe obtained geometries during the optimization process.
(a)
Iteration 2
(b)
Iteration 3
(c)
Iteration 5
(d)
Iteration 8
(e)
Iteration 11
(f)
Iteration 17


## 4.1.1. A Cavity with One Inlet and Two Outlets

We consider a cavity with one inlet and two outlets having the same section (see Figure2(a)) where Γ
in is a disc of center z
in
=
(
0,0.5,0.5
) and radius r
in
=
0.1, Γ
out
1 is a disc of center z
out
1
=
(
0.75,0.5,0
) and radius r
out
=
0.05, and Γ
out
2 is a disc of center z
out
2
=
(
0.75,0.5,1
) and radius r
out
=
0.05. The results of this case are given in Figures 3 and 4. The optimal geometry is presented in Figure 3(a). It is obtained for V
desired
=
0.15
|
Ω
|. A vertical cut of the obtained velocity is shown in Figure 3(b). The variation of 𝒥 is illustrated in Figure 3(c). We present in Figure 4 some geometries obtained during the optimization process.The initial domain and Figureelocity field.
(a)
The initial domain
(b)
The initial velocity fieldFirst case: cavity with one inlet and two outlets.
(a)
The optimal design
(b)
The obtained velocity field
(c)
Variation of the function𝒥The obtained geometries during the optimization process.
(a)
Iteration 2
(b)
Iteration 3
(c)
Iteration 6
(d)
Iteration 8
(e)
Iteration 13
(f)
Iteration 19


## 4.1.2. A Cavity with One Inlet and Three Outlets

In this case, we use a cavity with one inletΓ
in and three outlets Γ
out
1, Γ
out
2, and Γ
out
3 (see Figure 5(a)). The inlet and the outlets are defined as(i)
Γ
in is a disc of center (0,0.5,0.2) and diameter equal to 0.2;
(ii)
Γ
out
1 is a disc of center (1,0.5,0.4) and diameter equal to 0.066;
(iii)
Γ
out
2 is a disc of center (1,0.5,0.6) and diameter equal to 0.066;
(iv)
Γ
out
3 is a disc of center (1,0.5,0.8) and diameter equal to 0.066.The initial domain and velocity field.
(a)
The initial domain
(b)
The initial velocity fieldThe initial domain and the initial velocity field are given in Figure5. The results of this example are illustrated in Figures 6 and 7. The optimal domain is computed for V
desired
=
0.2
|
Ω
|. It is presented in Figure 6(a). A vertical cut (at y
=
0.5) of the velocity is shown in Figure 6(b). Figure 7 shows some geometries obtained during the optimization process.Second case: a cavity with one inlet and three outlets.
(a)
The optimal shape design
(b)
The obtained velocity fieldThe obtained geometries during the optimization process.
(a)
Iteration 2
(b)
Iteration 3
(c)
Iteration 5
(d)
Iteration 8
(e)
Iteration 11
(f)
Iteration 17


## 5. Conclusion

In this work, we have extended the topological gradient method for the nonstationary case. The discrete adjoint approach is introduced to overcome the difficulty coming from the nonlinearity of the operator. The proposed algorithm is applied to determine the optimal shape of tubes in a cavity. The optimal domain is obtained iteratively by inserting some obstacles in the initial one. The location and size of the obstacles are described by the topological gradient.The proposed approach has two main features. The first one concerns the adaptation for other nonstationary problems. The derived analysis is general and can be adapted for various operators like elasticity, Helmholtz, Maxwell, and so forthThe second interesting feature concerns the efficiency and the simplicity of the numerical algorithm. It is easy to be implemented and can be used for many applications. Only a few iterations are needed to construct the final domain. At each iteration, we only need to solve the direct and the adjoint problems on a fixed grid.



---

*Source: 923016-2014-03-24.xml*