Contents
Contents

1

8 Topology Optimization Problems of Density Variation Type
8.1 Set of Design Variables . . . . . . . . . . . . . . . . . . . . . . . .
8.2 State Determination Problem . . . . . . . . . . . . . . . . . . . .
8.3 Topology Optimization Problem of θ-Type . . . . . . . . . . . . .
8.4 Existence of an Optimum Solution . . . . . . . . . . . . . . . . .
8.5 Derivatives of Cost Functions . . . . . . . . . . . . . . . . . . . .
8.5.1 θ-Derivatives of Cost Functions . . . . . . . . . . . . . . .
8.5.2 Second-Order θ-Derivative of Cost Functions . . . . . . .
8.5.3 Second Order θ-Derivative of Cost Function Using
Lagrange Multiplier Method . . . . . . . . . . . . . . . . .
8.6 Descent Directions of Cost Functions . . . . . . . . . . . . . . . .
8.6.1 H 1 Gradient Method . . . . . . . . . . . . . . . . . . . . .
8.6.2 H 1 Newton Method . . . . . . . . . . . . . . . . . . . . .
8.7 Solution of Topology Optimization Problem of θ-Type . . . . . .
8.7.1 Gradient Method for Constrained Problems . . . . . . . .
8.7.2 Newton Method for Constrained Problems . . . . . . . . .
8.8 Error Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.9 Topology Optimization Problem of Linear Elastic Body . . . . .
8.9.1 State Determination Problem . . . . . . . . . . . . . . . .
8.9.2 Mean Compliance Minimization Problem . . . . . . . . .
8.9.3 θ-Derivatives of Cost Functions . . . . . . . . . . . . . . .
8.9.4 Second-Order θ-Derivatives of Cost Functions . . . . . . .
8.9.5 Second-Order θ-Derivative of Cost Function Using
Lagrange Multiplier Method . . . . . . . . . . . . . . . . .
8.9.6 Numerical Example . . . . . . . . . . . . . . . . . . . . .
8.10 Topology Optimization Problem of Stokes Flow Field . . . . . . .
8.10.1 State Determination Problem . . . . . . . . . . . . . . . .
8.10.2 Mean Flow Resistant Minimization Problem . . . . . . . .
8.10.3 θ-Derivatives of Cost Functions . . . . . . . . . . . . . . .
8.10.4 Second-Order θ-Derivatives of Cost Functions . . . . . . .
8.10.5 Second-Order θ-Derivative of Cost Function Using
Lagrange Multiplier Method . . . . . . . . . . . . . . . . .

3
7
8
11
13
17
17
19

1

22
24
25
28
29
29
31
32
38
38
39
40
41
44
45
49
49
51
52
54
56

2
8.10.6 Numerical Example . . . . . . . . . . . . . . . . . . . . .
8.11 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.12 Practice Problems . . . . . . . . . . . . . . . . . . . . . . . . . .

60
61
63

References

65

Index

69

Chapter 8Topology Optimization
Problems of Density
Variation Type
From this chapter, we ﬁnally examine shape optimization problems in continua.
Firstly, let us think about a problem seeking the appropriate arrangement ofholes in a domain where the boundary value problem of a partial diﬀerential
equation is deﬁned. Such a problem is known as the topology optimization
problem. Here, the term topology refers to the study of geometrical properties
and spatial relation of objects unaﬀected by the continuous change of their shapeor size. In mathematics, two mathematical objects are said to belong to thesame topology if they are images of two homotopic maps; that is, if one can becontinuously deformed into the other. Therefore, letting n be a natural number,a set of n-connected domains are regarded as belonging to the same homotopy
groups. Here, the term “topology optimization” in the topology optimization
problem refers to the determination of the connectivity of the design domain that
optimizes an object’s material distribution through insertion and arrangementof holes in its structure. However, as will be explained in detail later, the shapeof the holes actually becomes the target. Therefore, the problems dealt within this chapter also become included in shape optimization problems in a wider
sense. In this book, it will be referred to as the topology optimization problem
in the sense that topology is also in the scope of the design.
Until now, various problem formulations with respect to the topology
optimization problem and solutions to these problems have been proposed. By
setting up a ﬁxed domain in which a design target is included, a way to choose
the characteristic function (an L∞ -class function which takes the value 0 in holes
and 1 in a domain) of the domain as the design variable was considered. In sucha problem, we determine a state determination problem as a boundary value
problem multiplying the characteristic function with the coeﬃcient of the partial
diﬀerential equation, and deﬁne the cost function by the characteristic function
3

4

Chapter 8 Topology Optimization Problems of Density Variation Type

Fig. 8.1: Rank 2 material.
and the solutions of the boundary value problem. However, it was shown thatsuch a problem does not always have a solution [31]. The basic reason for this isthat a set of L∞ -class functions does not induce enough regularity to keep thecompactness of the admissible set of design variables [2]. Subsequently, the ideaof applying a method (homogenization technique) for describing a continuum
problem with a microstructure was shown [3, 19–22, 24, 25]. In particular, withrespect to a d ∈ {2, 3}-dimensional linear elastic body, the material constructedof microlayers of d types crossing over one another such as that in Fig. 8.1
are referred to as rank d material. Since the homogenized material constantof the rank d material could be analytically obtained, it was investigated ina lot of papers, see, for instance, [5, 7, 18, 32]. As a result, it was shown thatin a uniform stress ﬁeld, if each layer density in the direction of each principal
stress is determined in proportion to the value of each principal stress, the mean
compliance will be minimized as the volume is constrained [18]. However, noresults could be obtained which recognize macro holes. Moreover, since rank
d material has no rigidity with respect to shearing deformation, there was an
issue with it not being applicable as a mechanical structure.
The generation of holes was conﬁrmed when a continuum with a rectangularhole in a microcell Y similar to the one in Fig. 8.2 was assumed [6, 9, 26, 38].
This problem is constructed as a function optimization problem with the design⊤
variable as (a1 , a2 , θ) : D → R3 in Fig. 8.2. The numerical solution of thisproblem was obtained by an iterative method satisfying the optimality criteria[38]. If such a numerical solution is obtained, we can deﬁne the density by theratio of the magnitude of holes to the microcell and determine the shapes ofholes from the isosurface of the density with an appropriate threshold.
After that, it was shown that even if a microstructure is not assumed, whenfunction ϕ : D → [0, 1] (in reality, in order to avoid the discontinuity of thesolution of state determination problem occurring at ϕ → 0, ϕ : D → [c, 1] isused with some small constant c) representing density is set to be the design
variable, the same topology can be obtained as per the case of micro rectangularholes [28, 41]. In this case, a material characteristic (a coeﬃcient of a partial
diﬀerential equation) k was assumed to be given byk (ϕ) = k0 ϕα ,Chapter 8 Topology Optimization Problems of Density Variation Type

5

Fig. 8.2: A two-dimensional continuum with micro rectangular holes.

(a) Boundary value problem(b) Density ϕ and material constant k (ϕ)

Fig. 8.3: SIMP model.

where k0 is the existing material constant and α > 1 is a constant. Figure 8.3
shows the image of a topology optimization problem in that case. Figure 8.3(a) shows an example of a state determination problem (boundary value
problem). Figure 8.3 (b) shows the relationship between density ϕ and material
characteristic k. If the density which is the design variable ϕ (function deﬁnedon D) varies, the material characteristic k (ϕ) (function deﬁned on D) varies
via the function in Fig. 8.3 (b) and causes the solution of the boundary value
problem of Fig. 8.3 (a) to vary. A topology optimization problem deﬁned
in this way is referred to as topology optimization problem of density type.
Moreover, this problem is also referred to as the SIMP (solid isotropic material
with penalization) model [35]. The reason for this can be explained in thefollowing way. With respect to mid-level density (ϕ = 0.5), for materials withhomogeneous material such as in Fig. 8.4 (a), the material characteristic k
becomes 0.5α . But if, as in Fig. 8.4 (b), it splits into ϕ = 0 and ϕ = 1, k
becomes 0.5. For 0.5α < 0.5, it becomes a model which gives a penalty tomaterials split into ϕ = 0 and ϕ = 1 as having a greater material characteristic
value compared with uniform materials.
With respect to the topology optimization problem of density type, theFréchet derivatives of cost functions are calculated via the ﬁnite element method
using an evaluation method such as one shown later. However, if a constantdensity is assumed for each ﬁnite element and moved in the negative direction of
the Fréchet derivative, it has been pointed out that the density is seen to vibrate

6

Chapter 8 Topology Optimization Problems of Density Variation Type(a) Uniform material(b) Material split into 0 and 1

Fig. 8.4: Microstructures when density is 0.5.

(a) Checkerboard pattern(b) Optimal density by H 1 gradient method

Fig. 8.5: Numerical example of density-type topology optimization problem forlinear elastic body (provided by Quint Corporation).

in a checkerboard pattern [10]. Figure 8.5 (a) shows the results of numericalanalysis with respect to a mean compliance minimization problem (Problem
8.9.3) of a two-dimensional linear elastic body. In a state determination problem,
an external force pointing in the downward direction acts on the center of theright edge while the left edge is ﬁxed. The black and white elements show thatϕ is close to 1 and 0, respectively. When this type of phenomenon occurs, it can
be avoided by implementing post-processing such as ﬁltering, etc. [17, 23, 37].
Moreover, methods to approximate the distribution of material parameters in a
microstructure with continuous basis functions were shown [27]. However, it has
been pointed out that a diﬀerent issue called the island phenomenon, etc., mayarise [34]. Moreover, to remove the intermediate density regions, methods using
projection from intermediate density to zero one values were proposed [36, 40]．
In contrast, in this chapter, we shall think about the numerical analysis ofdensity-type topology optimization problem along the framework of abstract
optimization design problem shown in Chap. 7. However, in this chapter, afunction θ deﬁned on D is newly chosen as the design variable rather thanchoosing directly the density to be the design variable, for reasons shown later.
Hence, the density-type topology optimization problem in such a scenario will
be referred to as the topology optimization problem of θ-type. The frameworkof the logic used in this chapter is shown clearly in the paper [4]. Figure 8.5 (b)is the result obtained from the algorithm shown in Sect. 8.7. The fact that nonumerical instability phenomenon such as the checkerboard pattern is generated
can be seen.
This chapter is constructed in the following way. In Sect. 8.1, theadmissible set of the design variable θ is deﬁned in order to construct a topology

8.1 Set of Design Variables

7

optimization problem in continuum. In Sect. 8.2, a θ-type Poisson problem
is deﬁned as a state determination problem assuming that a design variable
is given. Using the design variable and the solution (state variables) of thestate determination problem, the topology optimization problem of θ-type will
be deﬁned in Sect. 8.3. Here, a cost function of general form will be used.
The existence of a solution to the topology optimization problem of θ-type isshown in Sect. 8.4. In Sect. 8.5, by referring to the Fréchet derivatives of costfunctions with respect to variation of the design variable θ as θ-derivatives, theprocess for obtaining the θ-derivatives and second-order θ-derivatives of costfunctions will be seen based on the methods for seeking Fréchet derivatives ofcost functions shown in Section 7.5. As a result, depending on the setting of
the state determination problem, it becomes apparent that the θ-derivatives of
the cost functions do not has regularity such that it would be included in theadmissible set of design variables. In Sect. 8.6, the abstract gradient method andabstract Newton method are speciﬁed with respect to the topology optimization
problem of θ-type. It becomes apparent that the variation of θ which can beobtained from such methods have the regularity such that it could be includedin the admissible set of design variable. In Sect. 8.7, the algorithm for solving
the topology optimization problem of θ-type will be considered. However, thebasic construction is the algorithm as shown in Section 3.7. Error estimations
when numerical analyses are conducted via this algorithm are considered in Sect.
8.8. Here, the results of error estimation from the numerical analysis shown in
Section 6.6 are used. Once the method for solving a Poisson-problem-related
topology optimization problem of θ-type has been conﬁrmed, we deﬁne a mean
compliance minimization problem of a linear elastic body and an energy loss
minimization problem of a Stokes ﬂow ﬁeld as topology optimization problems ofθ-type and show the process for seeking the θ-derivatives of their cost functionsin Sections 8.9 and 8.10. Moreover, a numerical example with respect to asimple problem is shown in each section.

8.1

Set of Design Variables

Firstly, let us deﬁne a set of design variables in order to construct a topology
optimization problem of a continuum. In this chapter, as shown in Fig. 8.3(a), D is taken to be a d ∈ {2, 3}-dimensional Lipschitz domain. ΓD ⊂ ∂D istaken to be a Dirichlet boundary and |ΓD | ̸= 0. ΓN ⊂ ∂D \ Γ̄D is a Neumannboundary.
In research so far, the range of density ϕ is limited to [0, 1]. The set offunctions with restricted range such as this cannot be a linear space. Hence inthis book, θ : D → R is set to be the design variable and the density is assumed
to be given by a sigmoid function ϕ ∈ C ∞ (R; R) with respect to θ. Several
functions are known to be sigmoid functions. Here, either
ϕ (θ) =

11
tan−1 θ +π2

(8.1.1)

8

Chapter 8 Topology Optimization Problems of Density Variation Typeorϕ (θ) =

1
1
tanh θ +2
2(8.1.2)is used. These graphs are shown in Fig. 8.6. At this point, there is a need tonote that ϕ : R → R is a function which returns (0, 1) when θ is given, and is
not a function deﬁned in D. However, because θ is a function with the domain
D, ϕ (θ) becomes a function deﬁned in D.
With respect to this type of design variable θ, let us use the framework ofan abstract optimal design problem (Problem 7.3.1) to deﬁne the linear spaceof design variables. As seen in Section 7.1, if the use of the gradient method is
considered, the linear space of design variables needs to be a real Hilbert space.
Hence, the linear space of design variable θ is set to be
(8.1.3)X = θ ∈ H 1 (D; R) θ = 0 in Ω̄C ,where Ω̄C ⊂ D̄ is a boundary or a domain on which a variation of θ is compressed
according to the design demands. If a function θC : D → R is speciﬁed and
it is assumed that θ = θC is established on Ω̄C , θ̃ = θ − θC is assumed to be
an element of X. In particular, if Ω̄C is not required, θ ∈ X = H 1 (D; R) is
assumed.
Furthermore, in order to be able to determine a Lipschitz continuous
boundary from the isosurfaces of θ and to be compact in X, we assume that theadmissible set of design variables is, at least, given byn
D = θ ∈ X ∩ H 2 (D; R) ∩ C 0,1 (D; R)oonmax ∥θ∥H 2 (D;R) , ∥θ∥C 0,1 (D;R) ≤ β(8.1.4)
where β is a positive constant. The requirement that D is a compact set in
X is assured by H 2 (D; R) ⋐ H 1 (D; R) obtained from the Rellich–Kondrachov
compact embedding theorem (Theorem 4.4.15). Moreover, in the same manneras Chap. 1, we consider that the boundedness constraint with norm is a side
constraint and assume that θ is an interior point of D (θ ∈ D◦ ) and when theside constraint is activated, we include it in the inequality constraints.

8.2

State Determination Problem

The linear space X and the admissible set D of design variables have been
deﬁned, hence let us next deﬁne the boundary value problem of a partial
diﬀerential equation which is a state determination problem. Here, the Poissonproblem is considered for simplicity.
A Poisson problem (Problem 5.1.1) is deﬁned in Chap. 5. Here, the Poissonproblem when θ is a design variable is called the θ-type Poisson problem andits deﬁnition is shown based on the framework of an abstract optimal design
problem (Problem 7.3.1).

8.2 State Determination Problem

9

1.0

0.5

0

2

4

6

(a) tan−1 θ/π + 1/2(b) (tanh θ + 1) /2

Fig. 8.6: Sigmoid functions for density ϕ (θ) with respect to the design variableθ.
The linear space of state variables (real Hilbert space) containing thehomogeneous solution (which is given by ũ = u − uD with a known function uD
provided for Dirichlet condition) of θ-type Poisson problem is set to be
(8.2.1)U = u ∈ H 1 (D; R) u = 0 on ΓD .
Furthermore, in order for the variation of θ obtained by the gradient method
which will be shown later to be included in D of Eq. (8.1.4), the admissible
set of state variables ũ of the homogeneous form with respect to the state
determination problem is set to beS = u ∈ U ∩ W 1,2qR (D; R) ∂ν u|ΓD ∈ L2 (ΓD ; R) ,(8.2.2)where we let qR be an integer satisfying qR > d.
In order for the homogeneous solution ũ with respect to a state determination
problem to be in S, from the results seen in Section 5.3, the following
assumptions are set with respect to the regularity of the known function.
Hypothesis 8.2.1 (Regularity of given functions) With respect to qR >d, it is assumed thatb ∈ C 1 X; L2qR (D; R) , pN ∈ W 1,2qR (D; R) , uD ∈ H 2 (D; R) ,where C 1 ( · ; · ) denotes the set of Fréchet diﬀerentiables (Deﬁnition 4.5.4). □Moreover, the following assumption is established with respect to theregularity of the boundary.
Hypothesis 8.2.2 (Opening angle of corner point)
Let D be a
two-dimensional domain. In relation to the corner points on the boundary, andwith respect to the Dirichlet boundary and Neumann boundary,(1) if the opening angle β is on the same type of boundary, β < 2π,(2) if it is on a mixed boundary, β < π.

10Chapter 8 Topology Optimization Problems of Density Variation Type

Fig. 8.7: θ-type Poisson problem.
Meanwhile, if D is a three-dimensional domain, the corner line on the boundaryis smooth and it is assumed that the aforementioned relationship holds onthe corner points of the boundary at a plane perpendicular to the corner line.
Furthermore, the crossing points of the corner lines or the apexes of the conicalboundaries are assumed not to have such singularities as they go beyond theframework within this book.
□
If Hypotheses 8.2.1 and 8.2.2 hold, the fact that u is included in
W 1,2qR (D; R) can be conﬁrmed as below. If Hypothesis 8.2.1 holds with respect
to given functions, as seen in Section 5.3.1, u is included in W 1,2qR class at all
points except those around the corner. Furthermore, from Proposition 5.3.1, if
ω >1−

2
2qR

(8.2.3)holds, u becomes included in the W 1,2qR class. Here, in the neighborhood aroundcorner points on the homogeneous boundary of Dirichlet or Neumann boundarytype, ω = π/β, hence the condition (1) in Hypothesis 8.2.2 can be obtained.
Moreover, at the neighborhood of a corner point of the mixed boundary, fromthe fact that ω = π/ (2β), the condition (2) of Hypothesis 8.2.2 is obtained.
Let us use the hypotheses above to deﬁne the state determination problem.
For simplicity, we consider a Poisson problem such as in Fig. 8.7. Here, ν · ∇ isalso to be written as ∂ν .
Problem 8.2.3 (θ-type Poisson problem) With respect to θ ∈ D, supposethat Hypotheses 8.2.1 and 8.2.2 are satisﬁed. Let α > 1 be a constant and ϕ (θ)
a function given by Eq. (8.1.1) or Eq. (8.1.2). In this case, obtain u : D → R
which satisﬁes−∇ · (ϕα (θ) ∇u) = b (θ)
α

in D,

ϕ (θ) ∂ν u = pN

on ΓN ,

u = uDon ΓD .
□The unique existence of a weak solution to Problem 8.2.3 is guaranteed by
the Lax–Milgram theorem (Theorem 5.2.4) with respect to ũ = u − uD ∈ U .

8.3 Topology Optimization Problem of θ-Type

11

Subsequently, ũ will be used to mean u − uD . Moreover, if Hypotheses 8.2.1 and8.2.2 are satisﬁed, u − uD ∈ S is guaranteed.
For later use, the Lagrange function with respect to Problem 8.2.3 is deﬁned
as
LS (θ, u, v)
Z=(−ϕα (θ) ∇u · ∇v + b (θ) v) dx
DZZ+pN v dγ +{(u − uD ) ϕα (θ) ∂ν v + vϕα (θ) ∂ν u} dγ, (8.2.4)
ΓN

ΓD

where u is not necessarily the solution of Problem 8.2.3, and v ∈ S was
introduced as a Lagrange multiplier. That v ∈ U is a Lagrange multiplierwith respect to Problem 8.2.3 can be conﬁrmed by recalling that in the process
to obtain the weak form of the Poisson problem (Problem 5.1.1) in Chap. 5,v ∈ U was introduced as a Lagrange multiplier. In Eq. (8.2.4), the third term on
the right-hand side is a term which was removed in Chap. 5 using u − uD , v ∈ U
when seeking the weak form of the Poisson problem. In reality, by removing
this term, u and v could be viewed as H 1 -class functions (in order for ∂ν v tohave meaning on ΓD , v needs to be a H 2 -class function). Here, however, thatterm will be left. The reason for this is because when a cost function fi includes
the boundary integral on ΓD , it becomes apparent that the boundary conditionof the adjoint problem with respect to fi is seen from the boundary integral on
ΓD in the Lagrange function of fi by using this term. Matching the deﬁnition
by Eq. (7.2.3) of a Lagrange function with respect to the abstract variational
problem in Chap. 7, using ũ = u − uD , we write
LS (θ, u, v) = −a (θ) (u, v) + l (θ) (v) = −a (θ) (ũ, v) + ˆl (θ) (v) ,where(8.2.5)Z
ϕα (θ) ∇u · ∇v dx,
ZDZ
pN v dγ,l (θ) (v) =b (θ) v dx +a (θ) (u, v) =D

(8.2.6)(8.2.7)ΓNˆl (θ) (v) = l (θ) (v) + a (θ) (uD , v) .

(8.2.8)When u is the solution to Problem 8.2.3,
LS (θ, u, v) = 0,
holds for all v ∈ U . This equation is equivalent to the weak form of Problem
8.2.3.

8.3

Topology Optimization Problem of θ-TypeThe design variable θ and solution u of the state determination problem (statevariable) were already deﬁned. Hence, let us use these to deﬁne the topology

12

Chapter 8 Topology Optimization Problems of Density Variation Type

optimization problem of θ-type. Here, we will consider a general cost function.
Let u be the solution of a state determination problem (Problem 8.2.3) withrespect to θ ∈ D and set the cost function for each i ∈ {0, 1, . . . , m} asZZfi (θ, u) =ζi (θ, u, ∇u) dx +
ηNi (u) dγDΓNZ−ηDi (ϕα (θ) ∂ν u) dγ − ci ,(8.3.1)
ΓD

where ci is a constant and is determined so that some (θ, ũ) ∈ D × S existsthat satisﬁes fi ≤ 0 for all i ∈ {1, . . . , m} (Slater constraint qualiﬁcation issatisﬁed). Moreover, ζi , ηNi and ηDi are assumed to be given by the following.
These hypotheses will be used in an adjoint problem (Problem 8.5.1) to satisfy
appropriate regularity requirements. To obtain the second-order θ derivativesof cost functions, additional hypotheses are needed but we will not specify
them further. Nevertheless, we will only assume that suﬃcient conditions are
satisﬁed by the state and adjoint state variables for us to be able to carry out
a second-order diﬀerentiation of the cost functions with respect to θ.
Hypothesis 8.3.1 (Regularity of cost functions)
For cost functionsfi(i ∈ {0, 1, . . . , m}) of Eq. (8.3.1), assume that ζi ∈ C 1 R × R × Rd ; R , ηNi ∈C 1 (R; R) and ηDi ∈ C 1 (R; R), and withrespect to (θ, u, ∇u, ∂ν u) ∈ D × S ×

G × GΓD (G = { ∇u | u ∈ S}, GΓD = ∂ν u|ΓD u ∈ S ),ζi (θ, u, ∇u) ∈ L1 (D; R) ,ζiθ (θ, u, ∇u) ∈ LqR (D; R) ,ζiu (θ, u, ∇u) ∈ L2qR (D; R) ,
ηNi (u) ∈ L1 (ΓN ; R) ,ζi(∇u)⊤ (θ, u, ∇u) ∈ W 1,2qR D; Rd ,′ηNi(u) ∈ L2 (ΓN ; R) ,

ηDi (ϕα (θ) ∂ν u) ∈ L1 (ΓD ; R) ,′ηDi(ϕα (θ) ∂ν u) ∈ W 1,2qR (ΓD ; R) .

□As a supplementary explanation, if we are worried that ∂ν u in the third termon the right-hand side of Eq. (8.3.1) cannot be deﬁned from the assumption ofu − uD ∈ S, it is remarked that this term will disappear due to the Dirichletcondition of adjoint problem (Problem 8.5.1) with respect to fi shown later.
Using cost functions f0 , f1 , . . . , fm of Eq. (8.3.1) and the framework ofan abstract optimal design problem (Problem 7.3.1), the topology optimization
problem of θ-type is deﬁned as follows.
Problem 8.3.2 (Topology optimization problem of θ-type) Let D and
S be given by Eq. (8.1.4) and Eq. (8.2.2), respectively. Suppose f0 , . . . ,fm : X × U → R are given by Eq. (8.3.1). In this case, obtain θ which satisﬁesmin(θ,u−uD )∈D×S{ f0 (θ, u) | f1 (θ, u) ≤ 0, . . . , fm (θ, u) ≤ 0, Problem 8.2.3} .
□8.3 Topology Optimization Problem of θ-Type

13

Based on the deﬁnition of a Lagrange function in Eq. (7.3.2) with respect
to the abstract optimal design problem in Chap. 7, the Lagrange function withrespect to Problem 8.3.2 is set to be
L (θ, u, v0 , v1 , . . . , vm , λ1 , . . . , λm )X
= L0 (θ, u, v0 ) +λi Li (θ, u, vi ) ,(8.3.2)i∈{1,...,m}⊤where λ = {λ1 , . . . , λm } ∈ Rm are Lagrange multipliers with respect to f1 ≤ 0,
. . . , fm ≤ 0. Moreover,Li (θ, u, vi ) = fi (θ, u) + LS (θ, u, vi )Z=
(ζi (θ, u, ∇u) − ϕα (θ) ∇u · ∇vi + b (θ) vi ) dx
DZ+(ηNi (u) + pN vi ) dγZΓN n+(u − uD ) ϕα (θ) ∂ν viΓDo+ (vi ϕα (θ) ∂ν u − ηDi (ϕα (θ) ∂ν u)) dγ − ci

(8.3.3)is a Lagrange function with respect to fi . Here, LS is the Lagrange function
with respect to Problem 8.2.3 deﬁned in Eq. (8.2.4). Moreover, vi is a Lagrange
multiplier with respect to a state determination problem prepared for fi and′assumes vi − ηDi∈ U . Furthermore, when thinking about the solution of the′
topology optimization problem of θ-type, the admissible set of ṽi = vi − ηDi(admissible set of adjoint variables) needs to be a subset of S.

8.4

Existence of an Optimum Solution

The existence of an optimum solution of Problem 8.3.2 can be assured by
Theorem 7.4.4 in Chap. 7. To use it, we need to show the compactness of
F = {(θ, ũ (θ)) ∈ D × S | Problem 8.2.3}(8.4.1)and the (lower semi) continuity of f0 . Here, we use ũ = u − uD ∈ U .
The compactness of F is presented in the following lemma.
Lemma 8.4.1 (Compactiness of F) Suppose that Hypothesis 8.2.1 andHypothesis 8.2.2 are satisﬁed. With respect to an arbitrary Cauchy sequenceθn → θ in X which is uniformly convergent in D and the solutions ũn = ũ (θn ) ∈U (n → ∞) of Problem 8.2.3, the convergenceũn → ũ strongly in U
holds, and ũ = ũ (θ) ∈ U solves Problem 8.2.3.

□

14Chapter 8 Topology Optimization Problems of Density Variation Type

Proof

Concerning the solution ũn of Problem 8.2.3 for θn , the inequality

αn ∥ũn ∥2U ≤ a (θn ) (ũn , ũn ) = ˆl (θn ) (ũn ) ≤ ˆl (θn )

U′∥ũn ∥U

holds, where a (θn ) and ˆl (θn ) are deﬁned in Eq. (8.2.5), and αn is a positive constantused in the deﬁnition of coerciveness for a (θn ) (see (1) in the answer to Exercise
5.2.5). When θn → θ (uniform convergence in D), αn can be replaced by a positive= ∥l (θn ) + a (θn ) (uD , · )∥ ′constant α not depending with n. The norm ˆl (θn )UU′

(l (θn ) is deﬁned in Eq. (8.2.5)) being bounded can be shown using (3) in the answer to
Exercise 5.2.5 via the convergence ϕα (θn ) → ϕα (θ) for θn → θ (uniform convergencel (v) and Ω in Exercise 5.2.5 are replaced by ˆl (θn ) (v) and D, respectively.
in D), where ˆHence, there exists a subsequence such that ũn → ũ weakly in U .
Next, we will show that ũ solves Problem 8.2.3 for θ. From the deﬁnition of
Problem 8.2.3, we havelim a (θn ) (ũn , v) = lim ˆl (θn ) (v) ,n→∞(8.4.2)n→∞with respect to an arbitrary v ∈ U . The right-hand side of Eq. (8.4.2) becomesl (θn ) (v) = ˆl (θ) (v) .
lim ˆ(8.4.3)n→∞Indeed, from Hypothesis 8.2.1, the inequalityZˆˆ(b (θn ) − b (θ)) v dx
l (θn ) (v) − l (θ) (v) =D≤ ∥b (θn ) − b (θ)∥L2 (D;R) ∥v∥L2 (D;R) → 0(n → ∞)

holds. The left-hand side of Eq. (8.4.2) becomes
lim a (θn ) (ũn , v) = a (θ) (ũ, v) .

(8.4.4)n→∞This is due to the fact that, since ũn → ũ weakly in U , we then have|a (θn ) (ũn , v) − a (θ) (ũ, v)|
ZZ=(ϕα (θn ) − ϕα (θ)) ∇ũn · ∇v dx +ϕα (θ) ∇ (ũn − ũ) · ∇v dx
DD≤ ∥ϕα (θn ) − ϕα (θ)∥C 0,1 (D;R) ∥ũn ∥H 1 (D;R) ∥v∥H 1 (D;R) + |a (θ) (ũn − ũ, v)|→0(n → ∞).

Substituting Eq. (8.4.3) and Eq. (8.4.4) into Eq. (8.4.2), the weak form of Problem
8.2.3 is obtained. It means that ũ = ũ (θ) ∈ U solves Problem 8.2.3.
Since the weak convergence of {un }n∈N to u was shown, the strong convergencecan be conﬁrmed by showing∥un ∥U → ∥u∥U(n → ∞).

Indeed, when we use a (θ) in Eq. (8.2.6) as a norm in U and deﬁne|||v||| = a (θ) (v, v) ,(8.4.5)

8.3 Topology Optimization Problem of θ-Type15

we have|||un ||| = a (θ) (un , un ) = a (θ − θn ) (un , un ) + a (θn ) (un , un )Z=
(ϕα (θ) − ϕα (θn )) ∇un · ∇un dx + l (θn ) (un )D→ l (θ) (u) = |||u|||(n → ∞).

(8.4.6)

From the above relation, it follows that un → u strongly in U , as desired.

□

We consider that the condition of ũ (θ) included in S is guaranteed in thesetting of Problem 8.2.3 (Hypotheses 8.2.1 and 8.2.2).
The continuity of f0 means that f0 is continuous on
(8.4.7)S = (θ, ũ (θ)) ∈ F | f1 (θ, u (θ)) ≤ 0, · · · , fm (θ, u (θ)) ≤ 0 .
Then, we will conﬁrm the continuity of f0 by showing the continuity of fi(i ∈ {0, 1, . . . , m}) by the following lemma and assuming that S is not empty.
Lemma 8.4.2 (Continuity of f0 ) Let fi be deﬁned as in Eq. (8.3.1) under
Hypothesis 8.3.1. Also, let un → u strongly in U which is determined by
Lemma 8.4.1 with respect to an arbitrary Cauchy sequence θn → θ in X, whichis uniformly convergent in D, satisfy ∥∂ν un − ∂ν u∥L2 (ΓD ;R) → 0 (n → ∞) on
ΓD . Then, fi is continuous with respect to θ ∈ D.
□
Proof

The proof will be completed whenZ(ζi (θn , un , ∇un ) − ζi (θ, u, ∇u)) dx|fi (θn , un ) − fi (θ, u)| ≤DZ(ηNi (un ) − ηNi (u)) dγ+ΓNZ(ηDi (ϕα (θn ) ∂ν un ) − ηDi (ϕα (θ) ∂ν u)) dγ+
ΓD= eD + eΓN + eΓD → 0(n → ∞)(8.4.8)is shown with respect to θn → θ (θn , θ ∈ D). From ũn → ũ weakly in
U , the convergence eD → 0 (n → ∞) is obtained. Indeed, writing ζ̃i (t) =ζi (tθn + (1 − t) θ, tu
 n + (1 − t) u, t∇un + (1 − t) ∇u) (t ∈ [0, 1]) and considering ζi ∈C 1 R × R × Rd ; R in Hypothesis 8.3.1, we get
ZZ

eD ≤ supζ̃iθ (t) [θn − θ] dx + supt∈[0,1]Zζ̃iu (t) [un − u] dxDζ̃i∇u (t) [∇un − ∇u] dx+ supt∈[0,1]≤ supt∈[0,1]DDζ̃iθ (t)t∈[0,1]+ supt∈[0,1]

LqR (D;R)ζ̃i∇u (t)∥θn − θ∥X + sup

W 1,2qR (D;Rd )t∈[0,1]ζ̃iu (t)L2qR (D;R)∥∇un − ∇u∥L2 (D;Rd ) ,∥un − u∥U

16

Chapter 8 Topology Optimization Problems of Density Variation TypeIn the same way, we can obtain eΓN → 0.
Meanwhile, the convergenceeΓD → 0 (n → ∞) can be shown in the following way. Writing η̃Di (t) =ηDi (ϕα (tθn + (1 − t) θ) (t∂ν un + (1 − t) ∂ν u)) (t ∈ [0, 1]) and considering ηDi ∈C 1 (R; R) in Hypothesis 8.3.1, we see that the sequence of inequalitiesZ′η̃Di(t) [ϕα (θn ) ∂ν un − ϕα (θ) ∂ν u] dxeΓD ≤ supt∈[0,1]

Z

ΓD
′
η̃Di
(t) [(ϕα (θn ) − ϕα (θ)) ∂ν un + ϕα (θ) (∂ν un − ∂ν u)] dx= supt∈[0,1]ΓD≤ ∥γΓD ∥3 supt∈[0,1]′η̃Di(t) W 1,2qR (D;R)

× ∥θn − θ∥X ∥∂ν un ∥L2 (ΓD ;R) + ∥θ∥X ∥∂ν un − ∂ν u∥L2 (ΓD ;R)holds. Here, boundedness of the trace operator ∥γΓD ∥ (Eq. (5.2.4)) was used. In
conclusion, we obtain Eq. (8.4.8).
□

From Lemma 8.4.1, the compactness of F was conﬁrmed. The continuityof f0 was shown by Lemma 8.4.2 and by the assumption that S is not empty.
Then, under these conditions, it can be assured that there exists an optimum
solution to Problem 8.3.2 by Theorem 7.4.4 (existence of an optimum solution).
Regarding the solution of Problem 8.3.2, we state the following remark.
Remark 8.4.3 (Existence of an optimum solution)
The compactnessof F deﬁned in Eq. (8.4.1) is based on the compactness of θ’sadmissible
set D deﬁned ino Eq. (8.1.4).
In Eq. (8.1.4), the conditionnmax ∥θ∥H 2 (D;R) , ∥θ∥C 0,1 (D;R) ≤ β with a positive constant β is added. Thiscondition corresponds to the condition called a side constraint in Chap. 1.
Such a side constraint is usually neglected, but it should be considered as a
non-equality constraint when the condition becomes active.
□
Moreover, regarding the selection of the linear space X and the admissible
set D for the design variable, the following remark is left for reference.
Remark 8.4.4 (Selection of X and D) The existence of a solution to
Problem 8.3.2 was conﬁrmed by Theorem 7.4.4 in Chap. 7. In the assumptionof the theorem, it is necessary to assume that D is a compact subset in X.
In this chapter, this relation was satisﬁed by taking X as a function spaceof class H 1 and D as a set of functions of (H 2 ∩ C 0,1 )-class based on the
Rellich–Kondrachov compact embedding theorem (Theorem 4.4.15). However,there are other selections. For instance, it is possible to select X as a function
space of C 0 class and D as a set of functions of C 0,1 class. In this case, theAscoli–Arzelà theorem (Theorem A.10.1) is used to show that D is a compact
subset in X [15, proof in Theorem 2.1, p. 16]. When those are selected, theassumptions and lemmas are changed to show the existence of a solution. In
this book, since a gradient method in a Hilbert space is considered, X and D
were selected as noted above.
□8.5 Derivatives of Cost Functions

8.5

17

Derivatives of Cost Functions

From this point onward, we will consider a solution to Problem 8.3.2 giventhat the conditions for its existence are satisﬁed. The Fréchet derivative of costfunction fi with respect to the variation of design variable θ will be referred toas a θ-derivative. Let us seek the θ-derivative of fi by the Lagrange multipliermethod such as that looked at in Section 7.5.2. Furthermore, let us seek thesecond-order θ-derivative of fi using a method such as that seen in Section 7.5.3.

8.5.1

θ-Derivatives of Cost Functions

Let us focus on the Lagrange function Li of fi deﬁned in Eq. (8.3.3). The
Fréchet derivative of Li with respect to an arbitrary variation (ϑ, û, v̂i ) ∈ X ×
U × U of (θ, u, vi ) becomesLi′ (θ, u, vi ) [ϑ, û, v̂i ]= Liθ (θ, u, vi ) [ϑ] + Liu (θ, u, vi ) [û] + Livi (θ, u, vi ) [v̂i ] .

(8.5.1)

The third term on the right-hand side of Eq. (8.5.1) becomes
Livi (θ, u, vi ) [v̂i ] = LSvi (θ, u, vi ) [v̂i ] = LS (θ, u, v̂i ) .

(8.5.2)Equation (8.5.2) is the Lagrange function of the state determination problem(Problem 8.2.3). Here, if u is the weak solution of the state determination
problem, the third term on the right-hand side of Eq. (8.5.1) is zero.
Moreover, the second term on the right-hand side of Eq. (8.5.1) becomes
Liu (θ, u, vi ) [û]Z −ϕα (θ) ∇vi · ∇û + ζiu û + ζi(∇u)⊤ · ∇û dx
=
D
Z
Z
′
′
+
ηNi û dγ +{ûϕα (θ) ∂ν v + (vi − ηDi) ϕα (θ) ∂ν û} dγ,ΓNΓD(8.5.3)′where
ζiu (θ, u, ∇u) [û],ζi(∇u)⊤ (θ, u, ∇u) [∇û],ηNi(u) [û]
and
α
′
′
α
′
ηDi (u) [ϕ (θ) ∂ν û] were written as ζiu û, ζi(∇u)⊤ · ∇û, ηNi û and ηDi ϕ (θ) ∂ν û,
respectively. Here, if vi is determined so that Eq. (8.5.3) becomes zero,
the second term on the right-hand side of Eq. (8.5.1) becomes zero. This
relationship is the weak form of the adjoint problem with respect to fi shown
next. Here, when vi is the weak solution of Problem 8.5.1, the second term of
the right-hand side of Eq. (8.5.1) vanishes. The boundary condition of Problem
8.5.1 is as shown in Fig. 8.8.

Problem 8.5.1 (Adjoint problem with respect to fi )For θ ∈ D◦ ,supposing the solution u is given to Problem 8.2.3, obtain vi : D → R which
satisﬁesin D,−∇ · (ϕα (θ) ∇vi ) = ζiu (θ, u, ∇u) − ∇ · ζi(∇u)⊤ (θ, u, ∇u)18

Chapter 8 Topology Optimization Problems of Density Variation Type

Fig. 8.8: Adjoint problem with respect to fi .
′ϕα (θ) ∂ν vi = ηNi(u) + ζi(∇u)⊤ · ν′
vi = ηDi

on ΓN ,

on ΓD .
□

Similarly to the solution u of the state determination problem, when′of Problem
Hypotheses 8.3.1 and 8.2.2 are satisﬁed, the solution ṽi = vi − ηDi
8.5.1 is guaranteed to be included in S.
Furthermore, the ﬁrst term on the right-hand side of Eq. (8.5.1) becomesLiθ (θ, u, vi ) [ϑ]Z=ζiθ + b′ vi − αϕα−1 ϕ′ ∇u · ∇vi ϑ dxDZ′+
αϕα−1 ϕ′ {(u − uD ) ∂ν vi + (vi − ηDi) ∂ν u} dγ.

(8.5.4)

ΓD

In the above equation, u and vi are assumed to be the weak solutions ofProblems 8.2.3 and 8.5.1, respectively. If we denote fi (θ, u) here by f˜i (θ), we
can writef˜i′ (θ) [ϑ] = Liθ (θ, u, vi ) [ϑ] = ⟨gi , ϑ⟩ ,(8.5.5)where
gi = ζiθ + b′ vi − αϕα−1 ϕ′ ∇u · ∇vi .

(8.5.6)

When Eq. (8.1.1) is used in ϕ (θ), we getϕ′ (θ) =1 1.
π 1 + θ2

(8.5.7)Moreover, when Eq. (8.1.2) is used,ϕ′ (θ) =1
1 1
1=
sech2 θ =2.
θ2
2 cosh2 θ(e + e−θ )(8.5.8)From the above, the following results can be obtained with respect toθ-derivative gi of fi .

8.5 Derivatives of Cost Functions

19

Theorem 8.5.2 (θ-derivative of fi ) For θ ∈ D◦ , suppose u and vi are theweak solutions of Problems 8.2.3 and 8.5.1 and these are said to be in S
of Eq. (8.2.2) (Hypotheses 8.2.1, 8.2.2 and 8.3.1 are satisﬁed). In this case,the θ-derivative of fi becomes Eq. (8.5.5). Hence, gi of Eq. (8.5.6) is in X ′ .
Furthermore, gi ∈ LqR (D; R).
□Proof
The fact that the θ-derivative of fi becomes the gi of Eq. (8.5.5) is as seen
above. The following results can be obtained in respect of the regularity of gi . If
Hölder’s inequality (Theorem A.9.1) and Poincaré’s inequality (Corollary A.9.4) areused in Eq. (8.5.5),|⟨gi , ϑ⟩|L1 (D;R)
≤ ∥ζiθ ∥LqR (D;R) + b′ L2qR (D;R) ∥vi ∥L2qR (D;R)

+ αϕα−1 ϕ′ C ∞ (R;R) ∥∇u∥L2qR (D;Rd ) ∥∇vi ∥L2qR (D;Rd ) ∥ϑ∥L2 (D;R)
≤ ∥ζiθ ∥LqR (D;R) + b′ L2qR (D;R) ∥vi ∥L2qR (D;R)
+ αϕα−1 ϕ′ C ∞ (R;R) ∥u∥W 1,2qR (D;R) ∥vi ∥W 1,2qR (D;R) ∥ϑ∥X

is obtained. The term (·) on the right-hand side of the equation above is completely
bounded due to the hypotheses. Hence, gi is included in X ′ . Moreover, from the factthat each term within (·) is in LqR (D; R), gi ∈ LqR (D; R) can be obtained.
□

From Theorem 8.5.2, the following can be said about the regularity of the
topology optimization problem of θ-type.
Remark 8.5.3 (Irregularity of topology optimization problem of θ-type)
If in Theorem 8.5.2, Hypotheses 8.2.1, 8.2.2 and 8.3.1 are made more strict anda problem is constructed such that u and vi are included in W 2,∞ (D; R), gi
would be included in C 0,1 (D; R). In this case, the design variable θ + ϵϑ (ϵ is apositive constant) updated via the gradient method such that −gi is replacedby ϑ would be included in the admissible set of design variables D. However,in such a case, it is necessary that the corner points permitted by Hypothesis
8.2.2 are removed, there is no boundary between the Dirichlet and Neumannboundaries, such as a Robin problem, or the neighborhoods of such points areincluded in Ω̄C in order to ﬁx θ.
If these conditions are not satisﬁed, gi is not included in C 0,1 (D; R). Hence,in a gradient method such that −gi is replaced by ϑ, θ + ϵϑ is not included in
the admissible set D of design variables. This result is thought to be one of thereasons for numerical instability phenomena in which a checkerboard pattern□
appears such as that in Fig. 8.5.

8.5.2

Second-Order θ-Derivative of Cost Functions

Furthermore, let us seek the second-order derivative (Hessian) of the costfunction with respect to the variation of the design variable. In Section 7.5.3,
the way to seek a second-order Fréchet derivative with respect to an abstract

20

Chapter 8 Topology Optimization Problems of Density Variation Type

optimal design problem has already been shown. Hence, let us follow thatmethod in order to seek the second-order θ-derivative of f˜i with respect to fi
given in Eq. (8.3.1).
The following assumption is established in order to obtain the second-orderθ-derivative of f˜i .
Hypothesis 8.5.4 (Second-order θ-derivative of f˜i ) With respect to thestate determination problem (Problem 8.2.3) and the cost function fi deﬁned
in Eq. (8.3.1), assume respectively that:(1) b is not a function of θ,(2) ζi is not a function of u (it is a function of θ and ∇u).
□Hypothesis 8.5.4 will be used in Eq. (8.5.16) to obtain Eq. (8.5.17). However,in the method shown in Sect. 8.5.3, this hypothesis will not be required.
The Lagrange function Li of fi is deﬁned by Eq. (8.3.3). Viewing (θ, u) as
a design variable and putting its admissible set and admissible set of directionsas
S = {(θ, u) ∈ D × S | LS (θ, u, v) = 0 for all v ∈ U } ,TS (θ, u) = {(ϑ, υ̂) ∈ X × U | LSθu (θ, u, v) [ϑ, υ̂] = 0 for all v ∈ U } ,the second-order Fréchet partial derivative of Li with respect to arbitrary
variations (ϑ1 , υ̂1 ) , (ϑ2 , υ̂2 ) ∈ TS (θ, u) of (θ, u) ∈ S, similarly to Eq. (7.5.21),
becomesLi(θ,u)(θ,u) (θ, u, vi ) [(ϑ1 , υ̂1 ) , (ϑ2 , υ̂2 )]= Liθθ (θ, u, vi ) [ϑ1 , ϑ2 ] + Liθu (θ, u, vi ) [ϑ1 , υ̂2 ]+ Liθu (θ, u, vi ) [ϑ2 , υ̂1 ] + Liuu (θ, u, vi ) [υ̂1 , υ̂2 ] .

(8.5.9)Each term on the right-hand side of Eq. (8.5.9) becomesZ′′Liθθ (θ, u, vi ) [ϑ1 , ϑ2 ] =ζiθθ − (ϕα (θ)) ∇u · ∇vi ϑ1 ϑ2 dx, (8.5.10)
DLiθu (θ, u, vi ) [ϑ1 , υ̂2 ]
Z n
o
′
=
ζiθ(∇u)⊤ · ∇υ̂2 − (ϕα (θ)) ∇υ̂2 · ∇vi ϑ1 dx,(8.5.11)D

Liuθ (θ, u, vi ) [ϑ2 , υ̂1 ]
Z n
o
′
=
ζiθ(∇u)⊤ · ∇υ̂1 − (ϕα (θ)) ∇υ̂1 · ∇vi ϑ2 dx,(8.5.12)D

Liuu (θ, u, vi ) [υ̂1 , υ̂2 ] = 0.

(8.5.13)′Here, the fact that u − uD , vi − ηDi, υ̂1 and υ̂2 become zero on ΓD was used.
Moreover,′(ϕα (θ)) = αϕα−1 (θ) ϕ′ (θ) ,(8.5.14)

8.5 Derivatives of Cost Functions

21′′(ϕα (θ)) = α (α − 1) ϕα−2 (θ) ϕ′2 (θ) + αϕα−1 (θ) ϕ′′ (θ) .

(8.5.15)On the other hand, with respect to arbitrary variations (ϑj , υ̂j ) ∈ TS (θ, u)
for j ∈ {1, 2}, the Fréchet partial derivative of LS becomesLSθu (θ, u, v) [ϑj , υ̂j ]Z′=− (ϕα (θ)) ϑ∇u − ϕα (θ) ∇υ̂j · ∇v dx
D(8.5.16)=0with respect to an arbitrary v ∈ U . Here, Hypothesis 8.5.4 and the fact that v
and υ̂j are both zero on ΓD were used. From Eq. (8.5.16), we can obtain∇υ̂j = −′(ϕα (θ))ϑj ∇u in D.
ϕα (θ)(8.5.17)

This relation becomes possible the following argument.
Substituting υ̂j of Eq. (8.5.17) into υ̂j in Eq. (8.5.9), and considering
Dirichlet boundary conditions in the state determination problem and the
adjoint problems with respect to f1 , . . . , fm , as well as υ̂j = 0 on ΓD , wehaveLiθu (θ, u, vi ) [ϑ1 , υ̂2 ] = Liuθ (θ, u, vi ) [υ̂1 , ϑ2 ]Z′o(ϕα (θ)) n α
′
(ϕ
(θ))
∇v
−
ζ
=
· ∇u ϑ1 ϑ2 dx.
⊤iiθ(∇u)αD ϕ (θ)(8.5.18)

Summarizing the results above, from Eq. (8.5.10), Eq. (8.5.13) andEq. (8.5.18), the second-order θ-derivative of f˜i becomeshi (θ, u, vi ) [ϑ1 , ϑ2 ]

Z 
′2
(ϕα (θ))′′α− (ϕ (θ)) ∇u · ∇vi
=
2 α
ϕ (θ)
D

′
(ϕα (θ))+ ζiθθ − 2 αζiθ(∇u)⊤ · ∇u ϑ1 ϑ2 dx
ϕ (θ)

Z 
ϕ′ (θ)=β (α, θ) ∇u · ∇vi + ζiθθ − 2αζiθ(∇u)⊤ · ∇u ϑ1 ϑ2 dx,
ϕ (θ)
D(8.5.19)where
β (α, θ) = α (α + 1) ϕα−2 (θ) ϕ′2 (θ) − αϕα−1 (θ) ϕ′′ (θ) .

(8.5.20)

When ϕ (θ) is given by Eq. (8.1.1) or Eq. (8.1.2), they respectively becomeβ (α, θ) = α (α + 1)

1
1
tan−1 θ +π2

α−2 1
π (1 + θ2 )

2

22

Chapter 8 Topology Optimization Problems of Density Variation Type(a) ϕ (θ) = tan−1 θ/π + 1/2(b) ϕ (θ) = (tanh θ + 1) /2

Fig. 8.9: Coeﬃcient function β (α, θ) in the second-order θ-derivative of the costfunction.
−α11
tan−1 θ +π2

α−1 (−)

2θ
2

π (1 + θ2 )(8.5.21)orα−2 
2
1
1
sech2 θtanh θ +
2
2
2

α−1

1
1
tanh θ +−α−sech2 θ tanh θ .
22β (α, θ) = α (α + 1)

(8.5.22)

Figure 8.9 shows the graph of β (α, θ). From these graphs, the fact thatβ (α, θ) > 0 holds can be conﬁrmed. Furthermore, if the remainder term in (·)
on the right-hand side of Eq. (8.5.19) is positive and bounded, hi (θ, u, vi ) [ · , · ]becomes a coercive and bounded bilinear form on X.

8.5.3

Second Order θ-Derivative of Cost Function Using
Lagrange Multiplier Method

When the Lagrange multiplier method is used to obtain the second-orderθ-derivative of a cost function, we use the same idea as proposed in Section 7.5.4.
Fixing ϑ1 , we deﬁne the Lagrange function with respect to f˜i′ (θ) [ϑ1 ] = ⟨gi , ϑ1 ⟩in Eq. (8.5.5) by
LIi (θ, u, vi , wi , zi ) = ⟨gi , ϑ1 ⟩ + LS (θ, u, wi ) + LAi (θ, vi , zi ) ,(8.5.23)

where LS is given by Eq. (8.2.4), andLAi (θ, vi , zi )Z =−ϕα (θ) ∇vi · ∇zi + ζiu zi + ζi(∇u)⊤ · ∇zi dx
D
Z
Z
′
′
+
ηNi zi dγ +{zi ϕα (θ) ∂ν v + (vi − ηDi) ϕα (θ) ∂ν zi } dγΓNΓD

(8.5.24)

8.5 Derivatives of Cost Functions

23

is the Lagrange function with respect to the adjoint problem (Problem 8.5.1)
with respect to fi . wi ∈ U and zi ∈ U are the adjoint variables provided for u
and vi in gi .
With respect to arbitrary variations (ϑ2 , û, v̂i , ŵi , ẑi ) ∈ X × U 4 of
(θ, u, vi , wi , zi ), the Fréchet derivative of LIi is written as
LIi′ (θ, u, vi , wi , zi ) [ϑ2 , û, v̂i , ŵi , ẑi ]= LIiθ (θ, u, vi , wi , zi ) [ϑ2 ] + LIiu (θ, u, vi , wi , zi ) [û]+ LIivi (θ, u, vi , wi , zi ) [v̂i ] + LIiwi (θ, u, vi , wi , zi ) [ŵi ]+ LIizi (θ, u, vi , wi , zi ) [ẑi ] .

(8.5.25)

The fourth term on the right-hand side of Eq. (8.5.25) vanishes if u is the solutionof the state determination problem. If vi can be determined as the solution of
the adjoint problem, the ﬁfth term of Eq. (8.5.25) also vanishes.
The second term on the right-hand side of Eq. (8.5.25) is
LIiu (θ, u, vi , wi , zi ) [û]
Z n

o

=
ζiθu + ζiuu zi + ζi(∇u)⊤ u · ∇zi û − αϕα−1 ϕ′ ∇vi · ∇û ϑ1D− ϕα ∇wi · ∇û dx.
(8.5.26)
Here, the condition that Eq. (8.5.26) is zero for arbitrary û ∈ U is equivalent tosetting wi to be the solution of the following adjoint problem.
Problem 8.5.5 (Adjoint problem of wi with respect to ⟨gi , ϑ1 ⟩) Under
the assumption of Problem 8.3.2, letting ϑ1 ∈ X be given, ﬁnd wi = wi (ϑ1 ) ∈ U
satisfying−∇ · (ϕα ∇wi ) = ∇ · αϕα−1 ϕ′ ∇vi + ζiθu + ζiuu zi + ζi(∇u)⊤ u · ∇zi ϑ1in D,
ϕ ∂ν wi = αϕα−1 ϕ′ ∂ν vi ϑ1αwi = 0on ΓN ,

on ΓD .
□The third term on the right-hand side of Eq. (8.5.25) is
LIivi (θ, u, vi , wi , zi ) [v̂i ]
Z
 ′

=
b v̂i − αϕα−1 ϕ′ ∇u · ∇v̂i ϑ1 − ϕα ∇zi · ∇v̂i dx.

(8.5.27)D

Here, the condition that Eq. (8.5.27) is zero for arbitrary v̂i ∈ U is equivalent
to setting zi to be the solution of the following adjoint problem.

24

Chapter 8 Topology Optimization Problems of Density Variation Type

Problem 8.5.6 (Adjoint problem of zi with respect to ⟨gi , ϑ1 ⟩) Under
the assumption of Problem 8.3.2, letting ϑ1 ∈ X be given, ﬁnd zi = zi (ϑ1 ) ∈ U
satisfying−∇ · (ϕα ∇zi ) = ∇ · αϕα−1 ϕ′ ∇u + b′ ϑ1 in D,
ϕα ∂ν zi = αϕα−1 ϕ′ ∂ν uϑ1
zi = 0on ΓN ,

on ΓD .
□

Finally, the ﬁrst term on the right-hand side of Eq. (8.5.25) becomesLIiθ (θ, u, vi , wi , zi ) [ϑ2 ]
Z h


=
ζiθθ u + b′′ vi − α (α − 1) ϕα−2 ϕ′2 + αϕα−1 ϕ′′ ∇u · ∇vi ϑ1Di− αϕα−1 ϕ′ (∇u · ∇wi + ∇vi · ∇zi ) + b′ (wi + zi ) ϑ2 dx.
Here, u, vi , wi (ϑ1 ) and zi (ϑ1 ) are assumed to be the weak solutions ofProblems 8.2.3, 8.5.1, 8.5.5 and 8.5.6, respectively. If we denote fi (θ, u) hereby f˜i (θ), we have the relation:
LIiθ (θ, u, vi , wi (ϑ1 ) , zi (ϑ1 )) [ϑ2 ] = f˜i′′ (θ) [ϑ1 , ϑ2 ]= ⟨gHi (θ, ϑ1 ) , ϑ2 ⟩ ,(8.5.28)

where the Hesse gradient gHi of fi is given bygHi (θ, ϑ1 )= − α (α − 1) ϕα−2 ϕ′2 + αϕα−1 ϕ′′ ∇u · ∇vi + ζiθθ u + b′′ vi ϑ1− αϕα−1 ϕ′ (∇u · ∇wi (ϑ1 ) + ∇vi · ∇zi (ϑ1 ))+ b′ (wi (ϑ1 ) + zi (ϑ1 )) .

(8.5.29)

We obtained two diﬀerent expressions for the second-order derivativef˜i′′ (θ) [ϑ1 , ϑ2 ]. Under the assumption of Hypothesis 8.5.4, they accord whenusing the same ϑ1 and ϑ2 . This relation will be conﬁrmed in Sections 8.9 and
8.10.

8.6

Descent Directions of Cost Functions

In Remark 8.5.3, it was shown that the topology optimization problem of θ-type
becomes irregular unless a special assumption of regularity is set. Here, let
us consider the gradient method and Newton method on the linear space X
of design variables with the functionality of regularizing θ-derivative of a costfunction. Here, with respect to the i ∈ {0, . . . , m}th cost function fi , assume
that the gradient gi ∈ X ′ of Eq. (8.5.6) and Hessian hi ∈ L2 (X × X; R) of
Eq. (8.5.19) are given and think about the method to obtain a descent directionof fi using the gradient method and Newton method on the linear space X ofdesign variables.

8.6 Descent Directions of Cost Functions

8.6.1

25

H 1 Gradient Method

The method for obtaining the aforementioned descent direction vector using
the solution ϑgi ∈ X to the next problem will be referred to as the H 1 gradientmethod of θ-type.
Problem 8.6.1 (H 1 gradient method of θ-type) Let X and D be
Eq. (8.1.3) and Eq. (8.1.4), respectively. Let aX : X × X → R be abounded and coercive bilinear form on X. In other words, it is supposed thatsome positive constants αX and βX exist and that2aX (ϑ, ϑ) ≥ αX ∥ϑ∥X ,|aX (ϑ, ψ)| ≤ βX ∥ϑ∥X ∥ψ∥X

(8.6.1)

holds with respect to arbitrary ϑ ∈ X and ψ ∈ X. For each fi ∈ C 1 (D; R), letgi (θk ) ∈ X ′ be its corresponding θ-derivative at θk ∈ D◦ which is not a local
minimum point. In this case, obtain ϑgi ∈ X which satisﬁesaX (ϑgi , ψ) = − ⟨gi (θk ) , ψ⟩with respect to an arbitrary ψ ∈ X.

(8.6.2)□

There is an arbitrary property for choosing aX : X × X → R which was
assumed in Problem 8.6.1. Several speciﬁc examples will be shown below.
Method Using the Inner Product of H 1 SpaceAn inner product in a real Hilbert space is coercive. Hence, let us use the innerproduct asZaX (ϑ, ψ) =(∇ϑ · ∇ψ + cD ϑψ) dx,(8.6.3)D

Here, it is assumed that cD is a uniformly bounded element of L∞ (D; R) whichis positive almost everywhere. In this case aX is a coercive bilinear form on X
(see solution to Exercise 5.2.7 (1)). Moreover, the following can be said for theway to choose cD . If cD takes a large value, the second term in the integral of
the right-hand side of Eq. (8.6.3) is dominant compared to the ﬁrst term andsuppresses the smoothing functionality. Hence, it becomes a result closer towhen −gi is chosen to be the direct search vector. Here, the size of the search
vector (step size) is considered to be adjusted by the size of the positive constantca used in the algorithm in Section 7.7.1 (same as Section 3.7).
The strong form of H 1 gradient method when using Eq. (8.6.3) is as below.
Problem 8.6.2 (H 1 gradient method using H 1 inner product) Let θ ∈D◦ , and gi ∈ X ′ of Eq. (8.5.6) be given. Find ϑgi : D → R which satisﬁes−∆ϑgi + cD ϑgi = −gi∂ν ϑgi = 0

in D,

on ∂D.
□

26

Chapter 8 Topology Optimization Problems of Density Variation Type

Fig. 8.10: H 1 gradient method using inner product in H 1 space.
Figure 8.10 shows the image of Problem 8.6.2. This problem is a boundary
value problem of an elliptic partial diﬀerential equation when Ω in the extended
Poisson problem of Problem 5.1.3 has changed to D, and c∂D = 0 is set. Hence,
numerical solutions can be obtained via numerical analysis methods such as theﬁnite element method.
Method Using Boundary Conditions
Moreover, even if the Dirichlet condition or Robin condition with respect to θ
are used, the bilinear form aX : X × X → R can be made coercive.
Firstly, let us think about using the Dirichlet boundary condition. On
Eq. (8.1.3) where the linear space X of the design variable is deﬁned, Ω̄C ⊂ D̄
was deﬁned to be a boundary or domain in which θ is ﬁxed under the design
demand. Here, the measure of the boundary or domain for Ω̄C is assumed tohave a positive value. In this case,ZaX (ϑ, ψ) =∇ϑ · ∇ψ dx(8.6.4)D\Ω̄C

is a bounded and coercive bilinear form on X as seen in the solution for Exercise
5.2.5. The strong form equation of H 1 gradient method in this case is as follows.
Problem 8.6.3 (H 1 gradient method using Dirichlet condition) Let
gi ∈ X ′ of Eq. (8.5.6) be given with respect to θ ∈ D◦ . Obtain ϑgi : D \ Ω̄C → R
which satisﬁes−∆ϑgi = −gi∂ν ϑgi = 0
ϑgi = 0in D \ Ω̄C ,on ∂D \ Ω̄C ,in Ω̄C .
□Problem 8.6.3 is a problem replacing Ω and ΓD in the Poisson problem
of Problem 5.1.1 with D \ Ω̄C and ∂ΩC , respectively. Figure 8.11 (a) shows
its image. The numerical solution of this problem can also be obtained via
numerical analysis method such as the ﬁnite element method.

8.6 Descent Directions of Cost Functions(a) Dirichlet condition

27(b) Robin condition

Fig. 8.11: H 1 gradient method using boundary conditions.
Furthermore, if the Robin condition is used, even if Ω̄C = ∅ is assumed in
Eq. (8.1.3), the coerciveness of aX (ϑ, ψ) can be obtained. Some positive-valued
and uniformly bounded function c∂D ∈ L∞ (∂D; R) is chosen andZZaX (ϑ, ψ) =∇ϑ · ∇ψ dx +c∂D ϑψ dγ.
(8.6.5)
D∂D

The fact that this aX becomes a coercive bilinear form in X is shown in thesolution of Exercise 5.2.7 (2). In this case, the strong form is as follows.
Problem 8.6.4 (H 1 gradient method using the Robin condition) Let
gi ∈ X ′ of Eq. (8.5.6) be given at θ ∈ D◦ . Obtain ϑgi : D → R which satisﬁes−∆ϑgi = −gi∂ν ϑgi + c∂D ϑgi = 0

in D,

on ∂D.
□

Figure 8.11 (b) shows an image of Problem 8.6.4. The Robin condition forthis problem is a condition typically used for a heat transfer boundary when aPoisson problem is viewed as a stationary heat transfer problem. The external
temperature at the boundary is set to zero and the heat transfer coeﬃcient is
set to be c∂D . Here, the numerical solutions of this problem can be obtained
via a numerical analysis method, such as the ﬁnite element method.
Regularity of H 1 Gradient Method
The following results can be obtained with respect to the weak solutions of
the H 1 gradient method (Problems 8.6.2 to 8.6.4) with respect to the topology
optimization problem of θ-type. Here, in this section, vicinities of the singularpoints as follows are denoted as B: when D is a two-dimensional domain,concave corner points on ∂D and corner points on ∂ΓD in mixed boundary
conditions for which the opening angle is greater than π/2, and when D isa three-dimensional domain, concave edges on ∂D and edges on ∂ΓD in mixed
boundary conditions for which the opening angle is greater than π/2. Moreover,fi (θ, u) when u is the solution to Problem 8.2.3 is written as f˜i (θ).

28

Chapter 8 Topology Optimization Problems of Density Variation Type

Theorem 8.6.5 (H 1 gradient method of θ-type) With respect to gi ∈LqR (D; R) in Theorem 8.5.2, the weak solutions ϑgi of Problems 8.6.2 to 8.6.4
exist uniquely, and ϑgi is in the H 2 (D; R) ∩ C 0,1 class on D \ B̄. Moreover, ϑgi
is a descent direction for this function.
□
Proof
From the fact that gi is in LqR (D; R) ⊂ X ′ , the Lax–Milgram theorem saysthat the weak solutions ϑgi of Problems 8.6.2 to 8.6.4 uniquely exist. Moreover, thefollowing results can be obtained regarding the regularity of the solution ϑgi . From
the fact that ϑgi satisﬁes an elliptic partial diﬀerential equation, the diﬀerentiability
increases by two orders compared to gi ; it becomes W 2,qR ⊂ H 2 (D; R) class on D \ B̄.
If Sobolev’s embedding theorem (Theorem 4.3.14) is applied to this, when qR > d,2−

d=1+σ >1qR

holds, where σ ∈ (0, 1). Therefore, in Theorem 4.3.14 (3), when p = qR , q = ∞, k = 1
and j = 1,
W 2,qR D \ B̄, R ⊂ C 0,1 D \ B̄, R
holds on D \ B̄. Then, ϑgi becomes H 2 (D; R) ∩ C 0,1 class on D \ B̄. Furthermore,with respect to the weak solutions ϑgi of Problems 8.6.2 to 8.6.4 and for some positive
constant ϵ̄, the estimatef˜i (θ + ϵ̄ϑgi ) − f˜i (θ) = ϵ̄ ⟨g i , ϑgi ⟩ + o (|ϵ̄|) = −ϵ̄aX (ϑgi , ϑgi ) + o (|ϵ̄|)≤ −ϵ̄αX ∥ϑgi ∥2X + o (|ϵ̄|)holds. Here, if ϵ̄ is taken to be suﬃciently small, f˜i (θ) decreases.

□

If the direction of variation of the design variable is determined using the H 1gradient method, a solution is found in the admissible set D of design variables
excepting the neighborhood of singular points from Theorem 8.6.5. From this,
it is thought that the H 1 gradient method is a regular gradient method.

8.6.2

H 1 Newton Method

Furthermore, if it is possible to calculate the second-order derivative (Hessian)hi ∈ L2 (X × X; R) of the cost function fi , the Newton method on X =H 1 (D; R) can be considered. Such a method is referred to as the H 1 Newton
method of θ-type.
Problem 8.6.6 (H 1 Newton method of θ-type)
Let X and D be
Eq. (8.1.3) and Eq. (8.1.4), respectively. For fi ∈ C 2 (D; R), its θ-derivativeand second-order θ-derivative at θk ∈ D◦ , which is not a local minimum point,
are taken respectively to be gi (θk ) ∈ X ′ and hi (θk ) ∈ L2 (X × X; R). Moreover,let aX : X ×X → R be a coercive and bounded bilinear form on X. Here, obtainϑgi ∈ X which satisﬁes
hi (θk ) [ϑgi , ψ] + aX (ϑgi , ψ) = − ⟨gi (θk ) , ψ⟩with respect to an arbitrary ψ ∈ X.

(8.6.6)□

8.7 Solution of Topology Optimization Problem of θ-Type29

In Problem 8.6.6, if the left-hand side of Eq. (8.6.6) is made to be just hi ,it cannot be expected to ﬁx the irregularity of gi (θk ) pointed out in Remark
8.5.3. In reality, hi calculated with Eq. (8.5.19) does not include the term of∇ϑ1 · ∇ϑ2 . Hence, in Problem 8.6.6, a bilinear form aX was added in order toensure coerciveness and boundedness of the left-hand side of Eq. (8.6.6) on X
and the regularity of ϑgi . For example, if Eq. (8.6.3) using an inner product on
X is to be used as a basis, let:ZaX (ϑ, ψ) =(cD1 ∇ϑ · ∇ψ + cD0 ϑψ) dx.
(8.6.7)D

Here, cD0 and cD1 are positive constants for achieving the coerciveness andregularity respectively. These have the same meaning as that explained after
Eq. (8.6.3).
Furthermore, in the case of the Newton method when the second-orderθ-derivative of fi (θ) is given by the Hesse gradient, Problem 8.6.6 is replacedwith the following problem.
Problem 8.6.7 (Newton method of θ-type using Hesse gradient) Let
X and D be Eq. (8.1.3) and Eq. (8.1.4), respectively. For fi ∈ C 2 (D; R),
the gradient of the θ-derivative of fi , search vector, which is obtained in theprevious step by the H 1 gradient method or H 1 Newton method of θ-typeusing the Hesse gradient, and Hesse gradient of fi at a non-local minimumpoint θk ∈ D◦ are denoted as gi (θk ) ∈ X ′ , ϑ̄gi ∈ X and gHi θk , ϑ̄gi ∈ X ′ ,respectively. aX : X × X → R is a coercive and bounded bilinear form on X.
Here, obtain a ϑgi ∈ X which satisﬁesaX (ϑgi , ψ) = − gi (θk ) + gHi θk , ϑ̄gi , ψ(8.6.8)with respect to an arbitrary ψ ∈ X.

8.7□

Solution of Topology Optimization Problemof θ-TypeThe abstract optimal design problem (Problem 7.3.1) and topology optimization
problem of θ-type (Problem 8.3.2) can be dealt with as in Table 8.1. Therefore
by appropriate replacements, the gradient method and Newton method withrespect to constrained problems shown in Section 7.7.1 (Section 3.7) and Section
7.7.2 (Section 3.8) can be applied.

8.7.1

Gradient Method for Constrained Problems

The gradient method with respect to a constrained problem can have a simple
algorithm such as Algorithm 3.7.2 shown in Section 3.7.1, which can be used byapplying changes such as those below:(1) Replace the design variable x and its variation y as θ and ϑ, respectively.

30

Chapter 8 Topology Optimization Problems of Density Variation Type

Table 8.1: Correspondence between abstract optimal design problem (Problem
7.3.1) and topology optimization problem of θ-type (Problem 8.3.2).

Design variable
State variableFréchet derivative of fi
Solution of gradient method

Abstract problem
ϕ∈Xu∈Ugi ∈ X ′φgi ∈ X

Topology optimization problemθ ∈ X = H 1 (D; R)u ∈ U = H 1 (D; R)
gi ∈ X ′ = H 1′ (D; R)ϑgi ∈ X = H 1 (D; R)(2) Equation (3.7.10) providing the gradient method is replaced by conditions
that establish
ca aX (ϑgi , ψ) = − ⟨gi , ψ⟩

(8.7.1)

with respect to an arbitrary ψ ∈ X, where aX (ϑgi , ψ) is a bilinear formon X used as the weak form of Problems 8.6.2 to 8.6.4.
(3) Replace Eq. (3.7.11) seeking the search vector withX
ϑg = ϑg0 +λi ϑgi ,nwhere IA =(8.7.2)i∈IAoi ∈ {1, . . . , m} | f˜i (θ) ≥ 0 .

(4) Replace Eq. (3.7.12) seeking the Lagrange multipliers with(⟨gi , ϑgj ⟩)(i,j)∈I 2 (λj )j∈IA = − (fi + ⟨gi , ϑg0 ⟩)i∈IA .

(8.7.3)A

Equation (8.7.3) is solved possibly several times, removing each time theconstraints where the associated Lagrange multiplier is negative (active
set method).
Furthermore, if a complicated algorithm such as Algorithm 3.7.6 is to be
used, the following changes should be added to (1) to (4) above:(5) The Armijo criterion Eq. (3.7.26) is replaced, with respect to ξ ∈ (0, 1),
with
+
*
X
L (θ + ϑg , λk+1 ) − L (θ, λ) ≤ ξ g0 +λi gi , ϑg .
(8.7.4)i∈IA(6) Replace the Wolfe criterion Eq. (3.7.27), with respect to µ (0 < ξ < µ < 1),
with
*
+
X
µ g0 +
λi gi , ϑgi∈IA

8.7 Solution of Topology Optimization Problem of θ-Type*≤

g0 (θ + ϑg ) +

X

31

+
λi k+1 gi (θ + ϑg ) , ϑg.

(8.7.5)i∈IA(7) Replace Eq. (3.7.21) for updating λk+1 based on the Newton–Raphson
method byδλ = (δλj )j∈IA−1= − (⟨gi (λk+1 l ) , ϑgj (λk+1 l )⟩)(i,j)∈I 2 (fi (λk+1 l ))i∈IA .

(8.7.6)A

Let us look at the points to be aware of when solving a topology optimization
problem of θ-type such as the one above.
In the topology optimization problem of θ-type (Problem 8.3.2), the solutionof the state determination problem with respect to a design variable θ ∈ X
or cost function becomes a non-convex non-linear mapping. This is because
the coeﬃcient ϕα (θ) of a partial diﬀerential equation used in a SIMP model
is a composite function of a sigmoid function and a power function. Hence,
depending on the deﬁnitions of cost functions and boundary conditions, there
may be cases when several local minimum points exist. In that case, the initial
distribution of θ needs to be changed and the convergence results need to becompared.
Moreover, the initial distribution of θ must be in the admissible set D
of design variables deﬁned in Eq. (8.1.4). In other words, it needs to be acontinuous function. If the initial distribution of θ is given by the characteristic
function (an L∞ class function) corresponding to the location of some holes,
caution needs to be taken that the discontinuities of θ at the boundaries of theholes are not removed, even when the methods above are used.
Furthermore, in topology optimization problems of θ-type shown in thischapter, a sigmoid function is used to change θ to ϕ. Therefore, as ϕ nears 0
and 1, there is a disadvantage that the gradient of ϕ with respect to θ becomes
small and convergence is slowed. This issue will hopefully to be improved using
the Newton method shown in the next section.

8.7.2

Newton Method for Constrained Problems

If the second-order θ-derivatives are computable in addition to the θ-derivativesof the cost functions, the gradient method with respect to a constrained problem
can be changed to a Newton method with respect to a constrained problem. In
this case, hi (θk ) [ϑgi , ψ] of Eq. (8.6.6) is replaced byX(8.7.7)hL (θk ) [ϑgi , ψ] = h0 (θk ) [ϑgi , ψ] +λik hi (θk ) [ϑgi , ψ] .
i∈IA (θk )In other words, let Eq. (8.6.6) bech hL (θk ) [ϑgi , ψ] + aX (ϑgi , ψ) = − ⟨gi (θk ) , ψ⟩ ,(8.7.8)

32

Chapter 8 Topology Optimization Problems of Density Variation Typewhere aX is deﬁned by Eq. (8.6.7), and ch , cD0 and cD1 are positive constants tocontrol the step size. Under this situation, a simple algorithm such as Algorithm
3.8.4 shown in Section 3.8.1 can be utilized after the following replacements:(1) Replace the design variable x and its variation y with θ and ϑ respectively.
(2) Replace Eq. (3.7.10) with the solution of Eq. (8.7.8).
(3) Replace Eq. (3.7.11) with Eq. (8.7.2).
(4) Replace Eq. (3.7.12) with Eq. (8.7.3).
When the second-order θ-derivative of fi (θ) is obtained as a Hesse gradient,Eq. (8.7.7) and Eq. (8.7.8) are replaced withXgHL θk , ϑ̄g = gH0 θk , ϑ̄g +λik gHi θk , ϑ̄g ,(8.7.9)i∈IA (θk )aX (ϑgi , ψ) = −gi (θk ) + ch gHL θk , ϑ̄g



,ψ ,

(8.7.10)respectively. Using the deﬁnitions, the following is added:(5) Replace Eq. (3.8.11) with Eq. (8.7.10).
Furthermore, when considering the complicated algorithm shown in Section
3.8.2, all sorts of innovations in response to the additional functionalities andcharacteristics of problems become necessary.

8.8

Error Estimation

If an algorithm such as the one shown in Sect. 8.7 is to be used to solve thetopology optimization problem of θ-type (Problem 8.3.2), the search vector ϑgcan be obtained via Eq. (8.7.2). In this regard, there is a need to obtain thesolutions to the boundary value problems of three elliptic partial diﬀerential
equations. In other words, the solution u to the state determination problem(Problem 8.2.3), the solutions v0 , vi1 , . . . , vi|I | of the adjoint problems
A(Problem 8.5.1) with respect to f0 , fi1 , . . . , fi|I | and the solutions ϑ0 , ϑi1 , . . . ,A

ϑi|I | from the H 1 gradient method of θ-type (Problem 8.6.1). Furthermore,A
there is a need to seek the Lagrange multipliers λi1 , . . . , λi|I | with Eq. (8.7.3).
A
Here, the numerical solutions with respect to the three boundary value problems
are assumed to be obtained by the ﬁnite element method, so let us use theresults of error estimation with respect to the numerical solutions from theﬁnite element method looked at in Section 6.6 in order to conduct the error
estimation of the search vector ϑg [29, 30].
Furthermore, if instead of the H 1 gradient method, the H 1 Newton method
is to be used, evaluations of the second-order derivatives of cost functions are
required. However, these will be omitted for now.

8.8 Error Estimation

33

In this section, D is assumed to be a polygon in two dimensions; a polyhedronin three dimensions and a regular ﬁnite element division T = {Di }i∈E withrespect to D is considered. Moreover, we deﬁne the maximum diameter h ofﬁnite elements as h (T ) of Eq. (6.6.2) and consider a sequence {Th }h→0 of ﬁnite
element divisions. Hereinafter, notation such as that below will be used:(1) Let the exact solutions of a state determination problem (Problem 8.2.3)
and adjoint problems (Problem 8.5.1) with respect to f0 , fi1 , . . . , fi|I | beAu and v0 , vi1 , . . . , vi|I | , respectively. Moreover, their numerical solutions
Afrom the ﬁnite element method can be written, with respect to i ∈ IA ∪{0},as
uh = u + δuh ,(8.8.1)

vih = vi + δvih .

(8.8.2)(2) Let the numerical solutions of θ-derivatives of cost functions f0 , fi1 , . . . ,fi|I | beA(8.8.3)gih = gi + δgih ,where gi and gih are functions of u, v0 , vi1 , . . . , vi|I | and uh , v0h , vi1 h ,A. . . , vi|I |h , respectively.
A(3) Let the exact solutions from the H 1 gradient method (for example,
Problem 8.6.2) calculated using the exact solutions g0 , gi1 , . . . , gi|I |Aof θ-derivatives be ϑg0 , ϑgi1 , . . . , ϑgi|I | . Moreover, the exact solutions
Afrom the H 1 gradient method calculated using the numerical solutions g0h ,gi1 h , . . . , gi|I | h are written, with respect to i ∈ IA ∪ {0}, asAϑ̂gi = ϑgi + δ ϑ̂gi .

(8.8.4)(4) Let the numerical solutions from the H 1 gradient method calculated using
the numerical solutions g0h , gi1 h , . . . , gi|I | h , with respect to i ∈ IA ∪ {0},Abeϑgih = ϑ̂gi + δ ϑ̂gih = ϑgi + δϑgih .

(8.8.5)(5) The coeﬃcient matrix (⟨gi , ϑgj ⟩)(i,j)∈I 2 of Eq. (8.7.3) constructed usingA
g0 , gi1 , . . . , gi|I | and ϑg0 , ϑgi1 , . . . , ϑgi|I | as A. Moreover, the coeﬃcientAAmatrix (⟨gih , ϑgjh ⟩)(i,j)∈I 2 of Eq. (8.7.3) constructed from g0h , gi1 h , . . . ,A
gi|I | h and ϑg0h , ϑgi1 h , . . . , ϑgi|I | h is denoted as Ah = A + δAh . Here,AA
we assume that fi = 0 and denote − (⟨gi , ϑg0 ⟩)i∈IA as b. Moreover,− (⟨gih , ϑg0h ⟩)i∈IA is denoted as bh = b + δbh . Furthermore, the exact
solution of the Lagrange multiplier is written as λ = A−1 b. Furthermore,
its numerical solution is written as
λh = (λih )i∈IA = A−1
h bh = λ + δλh .

(8.8.6)

34

Chapter 8 Topology Optimization Problems of Density Variation Type(6) Equation (8.7.2) constructed by λi1 h , . . . , λi|I | h and numerical solutionsAϑg0h , ϑgi1 h , . . . , ϑgi|I | h is written asAϑgh = ϑg0h +X

λih ϑgih = ϑg + δϑgh .

(8.8.7)i∈IA

In the deﬁnition above, δϑgh deﬁned in Eq. (8.8.7) represents the error of
the search vector. In this section, the aim is to evaluate the order of its norm∥δϑgh ∥X with respect to h. Here, the following hypothesis is established.
Hypothesis 8.8.1 (Error estimation of ϑg )In a state determination
problem and adjoint problems with respect to f0 , fi1 , . . . , fi|I | , let α > 1.
A
Moreover, the following is assumed with respect to qR > d and k, j ∈ {1, 2, . . .}:(1) The homogeneous form for the exact solutions u of a state determination
problem and v0 , vi1 , . . . , vi|I | of adjoint problems with respect to f0 , fi1 ,A. . . , fi|I | are elements of
A

Sk = U ∩ W k+1,2qR (D; R) .

(8.8.8)In order for this condition to be established, Hypotheses 8.2.1, 8.2.2 and8.3.1 need to be amended.
(2) The integrand of the cost function fi is, with respect to i ∈ IA ∪ {0},ζiθu ∈ L2qR (D; R) .
(3) There exist positive constants c1 , c2 , c3 which do not depend on h and fori ∈ IA ∪ {0},∥δuh ∥W j,2qR (D;R) ≤ c1 hk+1−j |u|W k+1,2qR (D;R) ,(8.8.9)∥δvih ∥W j,2qR (D;R) ≤ c2 h

(8.8.10)k+1−jδ ϑ̂gih

W j,2qR (D;R)|vi |W k+1,2qR (D;R) ,≤ c3 hk+1−j ϑ̂gi

W k+1,2qR (D;R),(8.8.11)is satisﬁed, where | · | expresses a semi-norm (see Eq. (4.3.12)).
(4) With respect to the coeﬃcient matrix Ah in Eq. (8.8.6), a positive constantc4 exists and
≤ c4
A−1
h
R|IA |×|IA |is satisﬁed, where ∥ · ∥R|IA |×|IA | represents the norm of a matrix (seeEq. (4.4.3)).
□

8.8 Error Estimation

35

In (1) of Hypothesis 8.8.1, since k ∈ {1, 2, . . .}, it is a stronger condition thanS deﬁned in Eq. (8.2.2). The reason for this is because in (3) of Hypothesis 8.8.1,
u and v0 , vi1 , . . . , vi|I | on the right-hand side of Eq. (8.8.9) and Eq. (8.8.10)
A

need to be in the W k+1,2qR class. Hypothesis 8.8.1 (3) is based on Corollary
6.6.4. Hypothesis 8.8.1 (4) is a condition which is established when gi1 , . . . ,gi|I | are linearly independent.
AHere, Theorem 8.8.5 shown later can be obtained. In order to show thisresult, the following three lemmas are used.
Lemma 8.8.2 (Error estimation of gi ) When the assumptions (1) and (2)in Hypothesis 8.8.1 as well as Eq. (8.8.9) and Eq. (8.8.10) are satisﬁed, withrespect to δgih of Eq. (8.8.3), there exists a positive constant c5 which does notdepend on h and the estimate⟨δgih , ϑ⟩ ≤ c5 hk ∥ϑ∥X
is established with respect to an arbitrary ϑ ∈ X.
Proof
δgih is yielded by the numerical error of δuh and δvih .
Eq. (8.5.5),□
Hence, from|⟨δgi , ϑ⟩| ≤ |Liθu (θ, u, vi ) [ϑ, δuh ] + Liθvi (θ, u, vi ) [ϑ, δvih ]|(8.8.12)is established. If with respect to the right-hand side of Eq. (8.8.12), Hölder’s inequality(Theorem A.9.1) and Poincaré’s inequality (Corollary A.9.4) are used,|Liθu (θ, u, vi ) [ϑ, δuh ] + Liθvi (θ, u, vi ) [ϑ, δvih ]|n≤ ∥ζiθu ∥L2qR (D;R) ∥δuh ∥L2qR (D;R) + b′ L2qR (D;R) ∥δvih ∥L2qR (D;R)+ αϕα−1 ϕ′ L∞ (D;R) ∥∇δuh ∥L2qR (D;Rd ) ∥∇vi ∥L2qR (D;Rd )o+ αϕα−1 ϕ′ L∞ (D;R) ∥∇uh ∥L2qR (D;Rd ) ∥∇δvih ∥L2qR (D;Rd ) ∥ϑ∥Xn≤ ∥ζiθu ∥L2qR (D;R) ∥δuh ∥W 1,2qR (D;R) + b′ L2qR (D;R) ∥δvih ∥W 1,2qR (D;R)+ αϕα−1 ϕ′ L∞ (D;R) ∥δuh ∥W 1,2qR (D;R) ∥vi ∥W 1,2qR (D;R)

o+ αϕα−1 ϕ′ L∞ (D;R) ∥uh ∥W 1,2qR (D;R) ∥δvih ∥W 1,2qR (D;R) ∥ϑ∥X

is established. Here, Eq. (8.8.9) and Eq. (8.8.10) in which j = 1 as well as Hypothesis
8.8.1 (1) leads to the result of the lemma.
□

Lemma 8.8.3 (Error estimation of ϑgi ) When Hypothesis 8.8.1 (1), (2)and (3) are satisﬁed, with respect to δϑgi of Eq. (8.8.5), the positive constantc6 which does not depend on h exists and the inequality∥δϑgih ∥X ≤ c6 hk
is established.

□

36

Chapter 8 Topology Optimization Problems of Density Variation TypeProof

The following is established from Eq. (8.8.4) and Eq. (8.8.5):∥δϑgih ∥X ≤ δ ϑ̂giX+ δ ϑ̂gih.

(8.8.13)X

represents the error of the exact solution of the H 1 gradient method (forHere δ ϑ̂gi
X

example, Problem 8.6.2) by δgih in Lemma 8.8.2 and δ ϑ̂gihwith respect to the numerical solution of H 1 gradient method.

represents the errorXδ ϑ̂giof Eq. (8.8.13)
X

satisﬁesaX δ ϑ̂gi , ϑ = − ⟨δgih , ϑ⟩with respect to an arbitrary ϑ ∈ X. Here, if ϑ = δ ϑ̂gi , the bound2

αX δ ϑ̂gi
X≤

DE
δgih , δ ϑ̂gi

(8.8.14)is established, where αX is a positive constant used in Eq. (8.6.1). If Lemma 8.8.2 isused with respect to δgih of Eq. (8.8.14), the estimateδ ϑ̂gi
X

≤

c5 k
h
αX

(8.8.15)

can be obtained. On the other hand, δ ϑ̂gihδ ϑ̂gihX≤ δ ϑ̂gih

W 1,2qR (D;R)

satisﬁes the inequalityX≤ c3 hk ϑ̂gi

(8.8.16)

W k+1,2qR (D;R)

from Eq. (8.8.11) in which j = 1. In Eq. (8.8.16),

ϑ̂gi

W k+1,2qR (D;R)

is bounded.

This is because if Hypothesis 8.8.1 (1) is used in the proof of Theorem 8.6.5, ϑ̂gi ∈W k+1,∞ (D; R) can be obtained. Hence, if Eq. (8.8.15) and Eq. (8.8.16) are substituted
into Eq. (8.8.13), the result for the lemma can be obtained.
□

Lemma 8.8.4 (Error estimation of λh ) When Hypothesis 8.8.1 is satisﬁed,
there exists a positive constant c7 which is not dependent on h, and the estimate∥δλh ∥R|IA | ≤ c7 hk
holds with respect to λh of Eq. (8.8.6).
Proof

□

With respect to λh of Eq. (8.8.6), the equation

δλh = A−1h (−δAh λ + δbh )
n 

= A−1
− (⟨δgih , ϑgj ⟩)(i,j)∈I 2 + (⟨gi , δϑgjh ⟩)(i,j)∈I 2 λhAAo+ (⟨δgih , ϑg0 ⟩)i∈IA + (⟨gi , δϑg0h ⟩)i∈IA
is established. If in Eq. (8.8.17), Hypothesis 8.8.1 (4) is used, the estimate∥δλh ∥R|IA |(8.8.17)

8.8 Error Estimation37≤ c4 1 + |IA | max |λi |
i∈IA

max

(i,j)∈IA ×(IA ∪{0})(|⟨δgih , ϑgj ⟩| + |⟨gi , δϑgjh ⟩|)(8.8.18)

holds. Since |IA | is bounded, with respect to |⟨δgih , ϑgj ⟩| of Eq. (8.8.18), the bound|⟨δgih , ϑgj ⟩| ≤ c5 hk ∥ϑgj ∥X(8.8.19)

is obtained from Lemma 8.8.2. Moreover, with respect to |⟨gi , δϑgjh ⟩|,|⟨gi , δϑgjh ⟩| ≤ c6 hk ∥gi ∥X(8.8.20)can be obtained from Lemma 8.8.3. In Eq. (8.8.20), ∥gi ∥X is bounded. This is because
if Hypothesis 8.8.1 (1) is used in the proof of Theorem 8.5.2, gi ∈ W k,qR (D; R) canbe obtained. Hence, if Eq. (8.8.18) and Eq. (8.8.19) are substituted into Eq. (8.8.17),
the result of the lemma can be obtained.
□

The following results can be obtained based on these lemmas.
Theorem 8.8.5 (Error estimation of ϑg )When Hypothesis 8.8.1 issatisﬁed, there exists a positive constant c not dependent on h, and∥δϑgh ∥X ≤ chk
is satisﬁed with respect to δϑgh of Eq. (8.8.7).
The following is established based on Eq. (8.8.7):Xδϑgh = δϑg0h +(δλih ϑgi + λi δϑgih ) .

□

Proof

(8.8.21)i∈IA

From Eq. (8.8.21),∥δϑgh ∥X ≤ 1 + |IA | max |λi |
i∈IA

max

i∈IA ∪{0}∥δϑgih ∥X+ ∥δλh ∥R|IA | |IA | max ∥ϑgi ∥Xi∈IA(8.8.22)can be obtained. If Eq. (8.8.22) is substituted with the results of Lemmas 8.8.3 and□8.8.4, the result of the theorem can be obtained.

From Theorem 8.8.5, the following can be said with respect to error
estimation of ﬁnite element solutions with respect to the topology optimization
problem of θ-type.
Remark 8.8.6 (Error estimation of ﬁnite element solution ϑgh ) Whenthe numerical solutions of the three boundary value problems (the state
determination problem, the adjoint problems and the H 1 gradient method) are
obtained by the ﬁnite element method with k = 1 order basis functions, from
Theorem 8.8.5, the error ∥δϑgh ∥X of search vector ϑgh reduces to the ﬁrst orderof h with respect to a sequence {Th }h→0 of ﬁnite element divisions.
□

38

Chapter 8 Topology Optimization Problems of Density Variation Type

Fig. 8.12: θ-type linear elastic body.

8.9

Topology Optimization Problem of Linear
Elastic Body

Let us change the state determination problem of θ-type topology optimization
problem to a linear elastic problem. Here, a mean compliance minimization
problem of a linear elastic body is deﬁned, and let us look at the θ-derivativeand second-order θ-derivative. If θ-derivatives and second-order θ-derivatives of
the cost functions can be obtained, such a problem can be solved in a similar
way to the Poisson problem.
Let D, ΓD and ΓN be the domain, Dirichlet boundary and Neumannboundary of a linear elastic problem, similar to the θ-type Poisson problem(Problem 8.2.3). For X and D, Eq. (8.1.3) and Eq. (8.1.4) respectively will be
used.

8.9.1

State Determination Problem

Let us deﬁne a linear elastic problem as a state determination problem. Let thelinear space U and admissible set S of state variables be
u = 0Rd on ΓD ,U = u ∈ H 1 D; Rd(8.9.1)
1,2qR
d
S =U ∩W
D; R .
(8.9.2)Moreover, Hypothesis 8.2.1 is changed in the following way.
Hypothesis 8.9.1 (Regularity of known functions) With respect to qR >d, assumeb ∈ C 1 X; L2qR D; Rd , pN ∈ L2qR ΓN ; Rd ,uD ∈ H 2 D; Rd , C ∈ L∞ D; Rd×d×d×d .
□On top of this, a problem such as the following is deﬁned with respect to aθ-type linear elastic body such as the one in Fig. 8.12.

8.9 Topology Optimization Problem of Linear Elastic Body

39

Problem 8.9.2 (θ-type linear elastic problem)
Let us suppose that
Hypotheses 8.9.1 and 8.2.2 hold. Moreover, let α > 1 be a constant and ϕ (θ)
is given by Eq. (8.1.1) or Eq. (8.1.2) with respect to θ ∈ D. In this case, obtainu : D → Rd satisfying−∇⊤ (ϕα (θ) S (u)) = b⊤ (θ)
α

in D,

(8.9.3)ϕ (θ) S (u) ν = pN

on ΓN ,

(8.9.4)u = uDon ΓD .

(8.9.5)□For later use, deﬁne the Lagrange function with respect to Problem 8.9.2 asZZαLS (θ, u, v) =(−ϕ (θ) S (u) · E (v) + b (θ) · v) dx +pN · v dγDΓNZ+{(u − uD ) · (ϕα (θ) S (v) ν) + v · (ϕα (θ) S (u) ν)} dγ,
ΓD

where u is not necessarily the solution of Problem 8.9.2. v ∈ U is a Lagrange
multiplier. If u is the solution of Problem 8.9.2,
LS (θ, u, v) = 0
holds with respect to an arbitrary v ∈ U .

8.9.2

Mean Compliance Minimization Problem

Let us deﬁne a topology optimization problem of θ-type with respect to a linear
elastic problem. Deﬁne the cost function as follows. With respect to the solution
u of Problem 8.9.2,
ZZZ
f0 (θ, u) =b (θ) · u dx +pN · u dγ −uD · (ϕα (θ) S (u) ν) dγ
D

ΓN

ΓD

(8.9.6)is referred to as the mean compliance. Moreover, the functionalZ
f1 (θ) =ϕ (θ) dx − c1

(8.9.7)

D

is referred to as the constraint function with respect to the domain measure of
the linear elastic body. Here, c1 is taken to be a positive constant, such that
f1 (θ) ≤ 0 holds with respect to some θ ∈ D. The reason that the functional f0
in Eq. (8.9.6) is called mean compliance is as follows. The ﬁrst and second termson the right-hand side of Eq. (8.9.6) are the work conducted by the volume force
b and traction pN , respectively. Since b and pN are ﬁxed, work conducted here
being small means that u is small. It can be referred to as external work if there
are only these two terms. However, the third term on the right-hand side of

40

Chapter 8 Topology Optimization Problems of Density Variation TypeEq. (8.9.6) is the negative value of the work done by uD . It is because the larger
work done by uD means that the resistance force with respect to deformationis stronger. Based on these, f0 will be called mean compliance in the sense that
it is the mean value for ease of deformation (compliance).
Using these deﬁnitions, the mean compliance minimization problem isdeﬁned as follows.
Problem 8.9.3 (Mean compliance minimization problem) Let D and S
be Eq. (8.1.4) and Eq. (8.9.2), respectively. Let f0 and f1 be Eq. (8.9.6) andEq. (8.9.7), respectively. In this case, obtain θ which satisﬁesmin(θ,u−uD )∈D×S{ f0 (θ, u) | f1 (θ) ≤ 0, Problem 8.9.2} .
□

8.9.3

θ-Derivatives of Cost Functions

Let us obtain the θ-derivative of f0 (θ, u) using the adjoint variable method. Let
the Lagrange function of f0 be
L0 (θ, u, v 0 )= f0 (θ, u) + LS (θ, u, v 0 )Z={−ϕα (θ) S (u) · E (v 0 ) + b (θ) · (u + v 0 )} dx
DZ+pN · (u + v 0 ) dγΓNZ+(u − uD ) · (ϕα (θ) S (v 0 ) ν)ΓD

+ (v 0 − uD ) · (ϕα (θ) S (u) ν) dγ.

(8.9.8)The Fréchet derivative of L0 with respect to arbitrary variation (ϑ, û, v̂ 0 ) ∈X × U × U of (θ, u, v 0 ) can be written as
L0′ (θ, u, v 0 ) [ϑ, û, v̂ 0 ]= L0θ (θ, u, v 0 ) [ϑ] + L0u (θ, u, v 0 ) [û] + L0v0 (θ, u, v 0 ) [v̂ 0 ] .

(8.9.9)Each term is considered below.
The third term on the right-hand side of Eq. (8.9.9) becomes
L0v0 (θ, u, v 0 ) [v̂ 0 ] = LSv0 (θ, u, v 0 ) [v̂ 0 ] = LS (θ, u, v̂ 0 ) .

(8.9.10)Equation (8.9.10) is a Lagrange function of the state determination problem(Problem 8.9.2). Hence, if u is the weak solution of the state determination
problem, the third term of the right-hand side of Eq. (8.9.9) is zero.
Moreover, the second term on the right-hand side of Eq. (8.9.9) becomes
L0u (θ, u, v 0 ) [û]

8.9 Topology Optimization Problem of Linear Elastic BodyZ

41

Z(−ϕ (θ) S (û) · E (v 0 ) + b (θ) · û) dx +pN · û dγ

α

=
D

ΓN

Z

{û · (ϕα (θ) S (v 0 ) ν) + (v 0 − uD ) · (ϕα (θ) S (û) ν)} dγ+
ΓD= LS (θ, v 0 , û) .

(8.9.11)Here, if v 0 is chosen so that Eq. (8.9.11) becomes zero, the second term on
the right-hand side of Eq. (8.9.9) vanishes. This relationship shows that the
self-adjoint relationshipu = v0(8.9.12)

holds.
Furthermore, the ﬁrst-term on the right-hand side of Eq. (8.9.9) becomesZ ′L0θ (θ, u, v 0 ) [ϑ] =b · (u + v 0 ) − αϕα−1 ϕ′ S (u) · E (v 0 ) ϑ dx.
D(8.9.13)Hence, u is taken to be a weak solution of Problem 8.9.2 and the self-adjoint
relationship (Eq. (8.9.12)) is assumed to hold. If f0 (θ, u) in this case is denotedas f˜0 (θ), we can writef˜0′ (θ) [ϑ] = L0θ (θ, u, v 0 ) [ϑ] = ⟨g0 , ϑ⟩ ,(8.9.14)where
g0 = 2b′ · u − αϕα−1 ϕ′ S (u) · E (u) .

(8.9.15)On the other hand, with respect to f1 (θ),
Z
f1′ (θ) [ϑ] =ϕ′ ϑ dx = ⟨g1 , ϑ⟩

(8.9.16)D

is established with respect to an arbitrary ϑ ∈ X.
Based on the above results, the function space which contains g0 of
Eq. (8.9.15) becomes the same result as Theorem 8.5.2. Hence, by applying H 1
gradient method, the fact that the search vector ϑg is in class C 0,1 is guaranteed.

8.9.4

Second-Order θ-Derivatives of Cost Functions

Furthermore, the second-order θ-derivatives of mean compliance f0 and theconstraint cost function f1 with respect to the domain measure of linear elastic
body can also be obtained. Here, we will follow the procedure shown in Sect.
8.5.2.
Firstly, let us think about the second-order θ-derivative of f0 . To correspondto Hypothesis 8.5.4 (1), here b is assumed not to be a function of θ. The
relationship corresponding to Hypothesis 8.5.4 (2) is satisﬁed here.

42

Chapter 8 Topology Optimization Problems of Density Variation TypeThe Lagrange function L0 of f0 is deﬁned by Eq. (8.9.8). Viewing (θ, u) as
a design variable, its admissible set and admissible set of directions are set as
S = {(θ, u) ∈ D × S | LS (θ, u, v) = 0 for all v ∈ U } ,TS (θ, u) = {(ϑ, υ̂) ∈ X × U | LSθu (θ, u, v) [ϑ, υ̂] = 0 for all v ∈ U } .
The second-order Fréchet partial derivative of L0 of Eq. (8.9.8) with respect toarbitrary variations (ϑ1 , υ̂ 1 ) , (ϑ2 , υ̂ 2 ) ∈ TS (θ, u) of design variable (θ, u) ∈ S
becomesL0(θ,u)(θ,u) (θ, u, v 0 ) [(ϑ1 , υ̂ 1 ) , (ϑ2 , υ̂ 2 )]= L0θθ (θ, u, v 0 ) [ϑ1 , ϑ2 ] + L0θu (θ, u, v 0 ) [ϑ1 , υ̂ 2 ]+ L0θu (θ, u, v 0 ) [ϑ2 , υ̂ 1 ] + L0uu (θ, u, v 0 ) [υ̂ 1 , υ̂ 2 ] .
Each term on the right-hand side of Eq. (8.9.17) becomesZ′′
L0θθ (θ, u, v 0 ) [ϑ1 , ϑ2 ] =− (ϕα (θ)) S (u) · E (v 0 ) ϑ1 ϑ2 dx,
DZ′
L0θu (θ, u, v 0 ) [ϑ1 , υ̂ 2 ] =− (ϕα (θ)) S (υ̂ 2 ) · E (v 0 ) ϑ1 dx,
DZ′
L0θu (θ, u, v 0 ) [ϑ2 , υ̂ 1 ] =− (ϕα (θ)) S (υ̂ 1 ) · E (v 0 ) ϑ2 dx,(8.9.17)(8.9.18)(8.9.19)(8.9.20)

D

L0uu (θ, u, v 0 ) [υ̂ 1 , υ̂ 2 ] = 0.

(8.9.21)

Here, u − uD , v 0 − uD , υ̂ 1 and υ̂ 2 use the fact that 0Rd on ΓD . Moreover,′′′(ϕα (θ)) and (ϕα (θ)) are Eq. (8.5.14) and Eq. (8.5.15), respectively. Here,with respect to an arbitrary variation (ϑj , υ̂ j ) ∈ TS (θ, u) for j ∈ {1, 2}, theFréchet partial derivative of Lagrange function LS of the state determination
problem becomesLSθu (θ, u, v) [ϑj , υ̂ j ]Z′=− (ϕα (θ)) ϑj S (υ) − ϕα (θ) S (υ̂ j ) · E (v) dx
D=0(8.9.22)with respect to an arbitrary v ∈ U . Here, the fact that v and υ̂ j are 0Rd on ΓD
is used. From Eq. (8.9.22), the equation
S (υ̂ j ) = −′(ϕα (θ))ϑj S (u)ϕα (θ)

in D(8.9.23)can be obtained. Hence, substituting υ̂ j of Eq. (8.9.23) into υ̂ 1 in Eq. (8.9.20)
and û2 in Eq. (8.9.19), and considering Dirichlet boundary conditions in thestate determination problem and the adjoint problems with respect to f1 , . . . ,fm as well as υ̂ j = 0Rd on ΓD , the equations
L0θu (θ, u, v 0 ) [ϑ1 , υ̂ 2 ] = L0θu (θ, u, v 0 ) [ϑ2 , υ̂ 1 ]8.9 Topology Optimization Problem of Linear Elastic Body(a) ϕ (θ) = tan−1 θ/π + 1/2

43(b) ϕ (θ) = (tanh θ + 1) /2

Fig. 8.13: Coeﬃcient functions ϕ′′ (θ) in h1 .
Z=′2(ϕα (θ))S (u) · E (v 0 ) ϑ1 ϑ2 dx
ϕα (θ)
D(8.9.24)

is obtained.
Summarizing the results above, by substituting Eq. (8.9.24) and Eq. (8.9.18)
into Eq. (8.9.17), the second-order θ-derivative of mean compliance f0 becomes)Z (′2(ϕα (θ))′′α2 α− (ϕ (θ)) S (u) · E (v 0 ) ϑ1 ϑ2 dx
h0 (θ, u, v 0 ) [ϑ1 , ϑ2 ] =
ϕ (θ)
D
Z
=
β (α, θ) S (u) · E (v 0 ) ϑ1 ϑ2 dx,(8.9.25)
D

where β (α, θ) is given by Eq. (8.5.20). Furthermore, if a self-adjoint relationshipis used, S (u)·E (v 0 ) > 0 and h0 (θ, u, v 0 ) [ · , · ] becomes a coercive and boundedbilinear form on X.
On the other hand, the second-order θ-derivative of f1 (θ) becomesZ′′
h1 (θ) [ϑ1 , ϑ2 ] = f1 (θ) [ϑ1 , ϑ2 ] =ϕ′′ (θ) ϑ1 ϑ2 dx(8.9.26)
Dwith respect to an arbitrary ϑ1 , ϑ2 ∈ X. Here, when ϕ (θ) of Eq. (8.1.1) is used,
we getϕ′′ (θ) = −12θ
.
π (1 + θ2 )2(8.9.27)

Moreover, when ϕ (θ) is given by Eq. (8.1.2), the equationϕ′′ (θ) = −sech2 θ tanh θ(8.9.28)

holds. Figure 8.13 shows the graphs of ϕ′′ (θ).
In this way, in the mean compliance minimization problem, the second-orderθ-derivative of the object cost function f0 is coercive but the second-orderθ-derivative of the constraint function f1 is not. Hence, if using the Newton
method (Problem 8.6.6), an additional term for capturing coerciveness becomes
necessary.

44

Chapter 8 Topology Optimization Problems of Density Variation Type

8.9.5

Second-Order θ-Derivative of Cost Function Using
Lagrange Multiplier Method

When the Lagrange multiplier method is used to obtain the second-orderθ-derivative of the mean compliance f0 , it becomes as follows. Fixing ϑ1 , wedeﬁne the Lagrange function for f˜0′ (θ) [ϑ1 ] = ⟨g0 , ϑ1 ⟩ in Eq. (8.9.14) by
LI0 (θ, u, w0 ) = ⟨g0 , ϑ1 ⟩ + LS (θ, u, w0 ) ,(8.9.29)

where LS is the Lagrange function of Problem 8.9.2, and w0 ∈ U is the adjoint
variable provided for u in g 0 .
With respect to arbitrary variations (ϑ2 , û, ŵ0 ) ∈ X × U 2 of (θ, u, w0 ), theFréchet derivative of LI0 is written as
LI0′ (θ, u, w0 ) [ϑ2 , û, ŵ0 ]= LI0θ (θ, u, w0 ) [ϑ2 ] + LI0u (θ, u, w0 ) [û]+ LI0w0 (θ, u, w0 ) [ŵ0 ] .

(8.9.30)

The third term on the right-hand side of Eq. (8.9.30) vanishes if u is the solutionof the state determination problem.
The second term on the right-hand side of Eq. (8.9.30) is
LI0u (θ, u, w0 ) [û]
Z
 ′

=
2b · û − 2αϕα−1 ϕ′ S (u) · E (û) ϑ1 − ϕα S (w0 ) · E (û) dx.
D(8.9.31)
Here, the condition that Eq. (8.9.31) is zero for arbitrary û ∈ U is equivalent
to setting w0 to be the solution of the following adjoint problem.
Problem 8.9.4 (Adjoint problem of w0 with respect to ⟨g0 , ϑ1 ⟩) Under
the assumption of Problem 8.9.2, let ϑ1 ∈ X be given. Find w0 = w0 (ϑ1 ) ∈ U
satisfying−∇⊤ (ϕα S (w0 )) = 2 ∇⊤ αϕα−1 ϕ′ S (u) + b′⊤ ϑ1 in D,
ϕα S (w0 ) ν = 2αϕα−1 ϕ′ S (u) νϑ1
w0 = 0Rd

on ΓN ,

on ΓD .
□

Finally, the ﬁrst term on the right-hand side of Eq. (8.9.30) becomes
LI0θ (θ, u, w0 ) [ϑ2 ]
Z h
 ′′

=
2b · u − α (α − 1) ϕα−2 ϕ′2 + αϕα−1 ϕ′′ S (u) · E (u) ϑ1Di− αϕα−1 ϕ′ S (u) · E (w0 ) + b′ · w0 ϑ2 dx.

8.9 Topology Optimization Problem of Linear Elastic Body

45(a) Initial density and boundary condition(b) H 1 gradient method (k = 100)(c) H 1 Newton method (k = 100)(d) H 1 Newton method(Hesse gradient, k = 100)

Fig. 8.14: Numerical example of mean compliance minimization: density
Here, u and w0 (ϑ1 ) are assumed to be the weak solutions of Problems 8.9.2
and 8.9.4, respectively. If we denote fi (θ, u) by f˜i (θ), we have the relation:
LI0θ (θ, u, v 0 , w0 (ϑ1 ) , z 0 (ϑ1 )) [ϑ2 ] = f˜0′′ (θ) [ϑ1 , ϑ2 ]= ⟨gH0 (θ, ϑ1 ) , ϑ2 ⟩ ,(8.9.32)

where the Hesse gradient gH0 of the mean compliance is given by
gH0 (θ, ϑ1 )= − α (α − 1) ϕα−2 ϕ′2 + αϕα−1 ϕ′′ S (u) · E (u) + 2b′′ · u ϑ1
− αϕα−1 ϕ′ S (u) · E (w0 (ϑ1 )) + b′ · w0 (ϑ1 ) .

(8.9.33)

If b is not a function of θ, with respect to the solution w0 of Problem 8.9.4,
E (w0 (ϑ1 )) = −2

αϕ′E (v 0 ) ϑ1
ϕ(8.9.34)

holds. Substituting Eq. (8.9.34) into Eq. (8.9.33), it can be conﬁrmed thatEq. (8.9.32) accords with Eq. (8.9.25).

8.9.6

Numerical Example

The results of mean compliance minimization for a two-dimensional linear elastic
body with a boundary condition referred to as the coat-hanging problem isshown in Figs. 8.14 to 8.16. The initial density (θ = 0) and the boundary
condition for the linear elastic problem are shown in Fig. 8.14 (a). A domain inwhich the density is constrained (Ω̄C0 in Eq. (8.1.3)) was not set. The programis written using the programming language FreeFEM (https://freefem.org/)[16] using the ﬁnite element method. In the ﬁnite element analyses of the linear

Chapter 8 Topology Optimization Problems of Density Variation Type

Cost function

46

1.0

1

0.8

0.8

0.6

0.6

0.4

0.4
0.2

0.2
0

20

40

60

80

100

(a) Cost functions

0

20

40

60

80

(b) Cost functions (search distance)
0

0

0

20

40

60

80

100

(c) Gradient of f0 on the search path

0

0.08

0.06

0.06

0.04

0.04

0.02

0.02
0

20

40

60

80

100

(e) Hessian of f0 on the search path

40

60

80

(d) Gradient of f0 on the search path(search distance)

0.08

0.00

20

0.00

0

20

40

60

80

(f) Hessian of f0 on the search path(search distance)

Fig. 8.15: Numerical example of mean compliance minimization: cost functions
and gradients and Hessians of f0 on the search path (gL : H 1 gradient method,hL , gL : H 1 Newton method, gH0 , h1 , gL : H 1 Newton method using the Hesse
gradient).

8.9 Topology Optimization Problem of Linear Elastic Body

47

4.0
3.5
4.0

3.0

3.5
2.5

3.0
2.5

2.0

2
1

2.0
1.5

20

40

60

80

100

(a) Iteration history

1.5

1.5

2.0

2.5

3.0

3.5

4.0

(b) (k − 1)-th vs. k-th plot

Fig. 8.16: Numerical example of mean compliance minimization problem:
distance ∥θk − θ∗ ∥X from an approximate minimum point θ∗ (gL : the gradientmethod, hL , gL : the Newton method，gH0 , h1 , gL : the Newton method using
the Hesse gradient).
elastic problem and the H 1 gradient method or the H 1 Newton method, thesecond-order triangular elements were used. In the case using the H 1 Newton
method, the routine of the H 1 Newton method was started from kN = 10. The
results were changed with ca in Eq. (8.7.1), cD in Eq. (8.6.3), cD1 and cD0 in
Eq. (8.6.7), ch in Eq. (8.7.8) and the parameter (errelas) to control the error
level in the adaptive mesh. For details, we refer the readers to the programs
themselves.1
Figure 8.14 (b) to (d) show the densities obtained by the three methods(H 1 gradient method using gL = g0 + λ1 g1 , H 1 Newton method using hL =h0 + λ1 h1 and gL , and H 1 Newton method using gH0 , h1 , gL ). Figure 8.15 (a)shows the cost functions f0 /f0init and 1 + f1 /c1 normalized with f0 at the initial
density denoted by f0init and c1 set with the domain integral of the initial density(volume), respectively, with respect to the iteration number k. Figure 8.15 (b)Pk−1shows those values with respect to the distance i=0 ϑg(i) X on the search
path in X. The graphs of f0 ’s gradient (the gradient of the Lagrange function
L = L0 + λ1 f1 ) calculated as gL , ϑg(k) / ϑg(k) X are shown in Fig. 8.15(c) and (d) with respect to the iteration number and the search distance,
respectively. Moreover, Fig. 8.15 (e) and (f) show the graphs of f0 ’s second-order
2derivative hL ϑg(k) , ϑg(k) / ϑg(k) X (in the case of the Newton method using2
the Hesse gradient, gH0 , ϑg(k) + λ1 h1 ϑg(k) , ϑg(k) / ϑg(k) X ) with respect
to the iteration number and the search distance, respectively. In these notations,the norm of the i-th search vector is deﬁned by
1/2
Z 

.
∇ϑg(i) · ∇ϑg(i) + ϑ2g(i) dx
ϑg(i) X =(8.9.35)
D

1 See Electronic supplementary material.

48

Chapter 8 Topology Optimization Problems of Density Variation TypeThe computational times until k = 100 by PC were 6.897, 17.073 and 32.394
sec by the H 1 gradient method, the H 1 Newton method and the H 1 Newton
method using the Hesse gradient, respectively.
We explain the numerical results and give some considerations as follows.
The graphs in Fig. 8.15 (a) clearly show that the convergence speed with respect
to the iteration number k is faster when using the H 1 Newton method than when
applying the H 1 gradient method. However, when the H 1 Newton method
started, cD1 and cD0 in Eq. (8.6.7) were replaced with smaller values so as
to avoid any numerical instability during iterations. As a result, it can beconsidered that the convergence speed was increased by enlarging the step size.
In reality, the following phenomenon was observed. When we set ch to zero(H 1 gradient method), the computation fails at k = kN . However, when weput a larger value for cD0 in Eq. (8.6.7), we have a convergence similar to the
H 1 Newton method. Moreover, based on the fact that the three graphs plottedin Fig. 8.15 (b) coincide, we conclude that the search paths due from eachmethods are actually the same. Based on these observations, we infer that the
H 1 Newton method is superior to the ﬁrst-order method, in the sense that we
can take larger values for the step size.
Based from the results shown in Fig. 8.15 (d) and (f), we also draw someparticular observations around the minimum point as follows. Firstly, sincethe Hessian of f0 on the search path has a positive value, we deduce that the
point of convergence is a local minimum. Secondly, we noticed that both theﬁrst derivative and the Hessian eventually diminish to zero. This key ﬁnding isobserved because we assumed a sigmoid function for the density of the design
valuable θ, and obtained a small variation of the density around the minimum
point where the density converges to 0 or 1 that causes a small variation of thecost functions.
In addition, Fig. 8.16 (a) shows the graphs of the distance θ(k) − θ∗ X froman approximate minimum point θ∗ obtained by the three methods with respect
to the iteration number k. The approximate minimum point θ∗ was given by thenumerical solution of θ when the iteration time is taken larger than the given
value in the H 1 Newton method. From this ﬁgure, it can be conﬁrmed that theconvergence orders for the results by the H 1 Newton methods are more than theﬁrst order. However, Fig. 8.16 (b), plotting the k-th distance ∥θk − θ∗ ∥X withrespect to the (k − 1)-th distance (the gradient of the graph shows the orderof convergence as explained by using Eq. (3.8.13)), shows that the convergenceorder of the H 1 Newton method is less than the second order, while more thanthe ﬁrst order. This result is possibly due to the fact that the bilinear form
aX in X was added to the original Hessian in order to ensure coercivenessand boundedness of the left-hand side of Eq. (8.6.6). By this addition, the H 1
Newton method has a diﬀerent structure from the original Newton method. It is
still unclear, however, whether a solution with second-order convergence exists
for the problem analyzed in this section.

8.10 Topology Optimization Problem of Stokes Flow Field

49

Fig. 8.17: Stokes ﬂow ﬁeld of θ-type.

8.10

Topology Optimization Problem of Stokes
Flow Field

The fact that the topology optimization problem can be constructed even withrespect to a ﬂow ﬁeld is shown in the literature [1, 8, 11–14, 33]. Here, the mean
ﬂow resistance minimization problem of the one-dimensional branched Stokes
ﬂow ﬁeld mentioned in Section 1.3 is extended to a d ∈ {2, 3}-dimensional
topology optimization problem of θ-type. Here, θ-derivatives of cost functions
and second-order θ-derivatives are shown up as far as can be sought. In thissection, D is assumed to be a Stokes ﬂow ﬁeld and X and D are taken to bedeﬁned in Eq. (8.1.3) and Eq. (8.1.4), repectively.

8.10.1

State Determination Problem

Let us deﬁne a Stokes problem as a state determination problem. A Stokes
problem (Problem 5.5.1) was deﬁned in Section 5.5 but here a Stokes ﬂow ﬁeld ofθ-type such as the one in Fig. 8.17 is considered. In this regard, some deﬁnitions
will be added. With respect to U and S, Eq. (8.9.1) and Eq. (8.9.2) will be used
respectively, but ΓD = ∂D. Furthermore, with respect to qR > d, putZq dx = 0 ,P = q ∈ L2 (D; R)(8.10.1)
D

Q = P ∩ L2qR (D; R) .

(8.10.2)

In a topology optimization problem of a ﬂow ﬁeld, a ﬂow ﬁeld passing
through porous media (penetration ﬂow) is used. In an penetration ﬂow,
between the ﬂow speed u and pressure p, the Darcy law given byku = − ∇pµis assumed to hold. Here, k and µ are positive constants known as penetrationand viscosity coeﬃcients. In a topology optimization problem of a ﬂow ﬁeld,
replace the constant µ/k representing the diﬃculty of penetration withα (1 − ϕ)ϕ (1 + α)(8.10.3)= ψ1
ψ (ϕ) = ψ1 1 −
ϕ+α
α+ϕ

50

Chapter 8 Topology Optimization Problems of Density Variation Type

Fig. 8.18: Coeﬃcient ψ expressing the ﬂow resistance with respect to ﬂuid
content ϕ.
and ∇p in the Stokes equation with ψ (ϕ (θ)) u + ∇p. Here, ϕ represents theﬂuid content equivalent to the density of ﬂuid and its range is assumed to belimited to [0, 1]. Hence, similarly to the density in the previous sections, ϕ isassumed to be given by the sigmoid function with respect to the design variableθ ∈ X. Moreover, ψ1 is a positive constant which gives the maximum value ofthe resistance to ﬂow, α is a constant controlling the non-linearity and is chosenfrom (0, 1]. In the paper [1], it is changed from 0.01 to 1 along with calculation
progression. Figure 8.18 shows the function ψ (ϕ).
Here, the following assumption is set.
Hypothesis 8.10.1 (Regularities of known functions) With respect toqR > d,b ∈ L2qR D; Rd , uD ∈ u ∈ W 1,2qR D; Rd∇ · u = 0 in D□

is assumed.

Using these assumptions, a Stokes problem of θ-type is deﬁned in thefollowing way.
Problem 8.10.2 (θ-type Stokes problem) Let b and uD satisfy Hypothesis
8.10.1, and Hypothesis 8.2.2 holds with respect to the opening angles ofboundary corner points.
Moreover, ψ (ϕ) is taken to be Eq. (8.10.3).
Furthermore, ϕ (θ) is assumed to be given by Eq. (8.1.1) or Eq. (8.1.2). Here,obtain (u, p) : D → Rd+1 which satisﬁes−∇⊤ µ∇u⊤ + ψ (ϕ (θ)) u⊤ + ∇⊤ p = b⊤ in D,(8.10.4)∇·u=0Z

u = uDp dx = 0.

in D,
on ∂D,(8.10.5)(8.10.6)(8.10.7)

D□8.10 Topology Optimization Problem of Stokes Flow Field

51

For later use, the Lagrange function with respect to Problem 8.10.2 is deﬁned
as
LS (θ, u, p, v, q)
Z



=
−µ ∇u⊤ · ∇v ⊤ − ψ (ϕ (θ)) u · v
D+ p∇ · v + b · v + q∇ · u dx
Z+{(u − uD ) · (µ∂ν v − qν) + v · (µ∂ν u − pν)} dγ,(8.10.8)∂D

where (u, p) is not necessarily the solution of Problem 8.10.2, and (v, q) ∈ U ×P
is a Lagrange multiplier.
When (u, p) is a solution to Problem 8.10.2, the equation
LS (θ, u, p, v, q) = 0
is established with respect to an arbitrary (v, q) ∈ U × P .

8.10.2

Mean Flow Resistant Minimization Problem

Let us deﬁne a topology optimization problem of θ-type with respect to a Stokes
ﬂow ﬁeld. Deﬁne the cost functions as follows. Firstly, let us deﬁne a costfunction representing the resistance to ﬂow asZZ
f0 (θ, u, p) = −b · u dx +uD · (µ∂ν u − pν) dγ.
(8.10.9)
D∂D

The ﬁrst term on the right-hand side of Eq. (8.10.9) represents the negative valueof the rate of work due to the volume force. This value was given a negative sign
because the greater it is, the greater the ﬂow speed is. On the other hand, thesecond term on the right-hand side of Eq. (8.10.9) is equivalent to the energy
per unit time lost through the viscosity inside Stokes ﬂow ﬁeld represented by
the boundary integral. From the fact that these express the property of ﬂow
resistance, f0 will be referred to as the mean ﬂow resistance. With respect tothis,Z
f1 (θ) =ϕ (θ) dx − c1
(8.10.10)
Dis the cost function for constraint with respect to the domain measure of theﬂow ﬁeld. Here, c1 is a positive constant such that f1 (θ) ≤ 0 holds with respect
to some θ ∈ D. Using these, the minimization problem of mean ﬂow resistanceis deﬁned as follows.
Problem 8.10.3 (Mean ﬂow resistance minimization problem) Let D,
S and Q be Eq. (8.1.4), Eq. (8.9.2) and Eq. (8.10.2), respectively. Let (u, p) be

52

Chapter 8 Topology Optimization Problems of Density Variation Typethe solution of Problem 8.10.2 with respect to θ ∈ D and f0 and f1 are given
by Eq. (8.10.9) and Eq. (8.10.10). In this case, obtain θ which satisﬁesmin(θ,u−uD ,p)∈D×S×Q{ f0 (θ, u, p) | f1 (θ) ≤ 0, Problem 8.10.2} .
□8.10.3θ-Derivatives of Cost Functions

Let us obtain the θ-derivative of f0 (θ, u, p) via the adjoint variable method. Let
the Lagrange function of f0 be
L0 (θ, u, p, v 0 , q0 )= f0 (θ, u, p) − LS (θ, u, p, v 0 , q0 )
Z n


=
µ ∇u⊤ · ∇v ⊤0 + ψ (ϕ (θ)) u · v 0 − p∇ · v 0Do− b · (v 0 + u) − q0 ∇ · u dx
Z−{(u − uD ) · (µ∂ν v 0 − q0 ν) + (v 0 − uD ) · (µ∂ν u − pν)} dγ.
∂D(8.10.11)
Comparing with Eq. (8.9.8) deﬁning the Lagrange function with respect to mean
compliance of linear elastic body, here a negative sign was put on LS . Thischange is so that a self-adjoint relationship can be obtained later. Although
in the mean compliance minimization problem of a linear elastic body, theminimization of the displacement was the aim, in a Stokes ﬂow ﬁeld mean ﬂow
resistance minimization problem, the maximization of ﬂow is aimed for. Hence,
this type of diﬀerence arose. The Fréchet derivative of L0 with respect to an
2
arbitrary variation (ϑ, û, p̂, v̂ 0 , q̂0 ) ∈ X × (U × P ) of (θ, u, p, v 0 , q0 ) isL0′ (θ, u, p, v 0 , q0 ) [ϑ, û, p̂, v̂ 0 , q̂0 ]= L0θ (θ, u, p, v 0 , q0 ) [ϑ] + L0up (θ, u, p, v 0 , q0 ) [û, p̂]+ L0v0 q0 (θ, u, p, v 0 , q0 ) [v̂ 0 , q̂0 ] .

(8.10.12)

Each term is considered below.
The third term on the right-hand side of Eq. (8.10.12) becomes
L0v0 q0 (θ, u, p, v 0 , q0 ) [v̂ 0 , q̂0 ] = LSv0 q0 (θ, u, p, v 0 , q0 ) [v̂ 0 , q̂0 ]= −LS (θ, u, p, v̂ 0 , q̂0 ) .

(8.10.13)

Equation (8.10.13) is a Lagrange function of the state determination problem(Problem 8.10.2). Hence, if (u, p) is the weak solution of the state determination
problem, the third term on the right-hand side of Eq. (8.10.12) vanishes.
Moreover, the second term on the right-hand side of Eq. (8.10.12) becomes
L0up (θ, u, p, v 0 , q0 ) [û, p̂]8.10 Topology Optimization Problem of Stokes Flow Field

53

Z n 


=
µ ∇û⊤ · ∇v ⊤0 + ψ (ϕ (θ)) û · v 0 − p̂∇ · v 0Do− b · û − q0 ∇ · û dx
Z−{û · (µ∂ν v 0 − q0 ν) + (v 0 − uD ) · (µ∂ν û − p̂ν)} dγ∂D= −LS (θ, v 0 , q0 , û, p̂)(8.10.14)

with respect to an arbitrary variation (û, p̂) ∈ U × P of (u, p). Hence, when the
self-adjoint relationship(u, p) = (v 0 , q0 )(8.10.15)

holds, the second term on the right-hand side of Eq. (8.10.12) becomes zero.
Furthermore, the ﬁrst term on the right-hand side of Eq. (8.10.12) becomesZ
L0θ (θ, u, p, v 0 , q0 ) [ϑ] =ψ ′ (ϕ (θ)) ϕ′ (θ) u · v 0 ϑ dx.

(8.10.16)

D

Therefore, suppose (u, p) is a weak solution of Problem 8.10.2 and that the
self-adjoint relationship Eq. (8.10.15) holds. If f0 (θ, u, p) in this case is written
as f˜0 (θ), the equationf˜0′ (θ) [ϑ] = L0θ (θ, u, p, v 0 , q0 ) [ϑ] = ⟨g0 , ϑ⟩

(8.10.17)

holds, whereg0 = ψ ′ ϕ′ u · u.

(8.10.18)

When using ψ (ϕ) of Eq. (8.10.3), we getψ ′ (ϕ) = −ψ1

α (1 + α)

2 .

(ϕ + α)(8.10.19)

On the other hand, with respect to f1 (θ),
f1′ (θ) [ϑ] =Zϕ′ ϑ dx(8.10.20)

D

holds with respect to an arbitrary ϑ ∈ X.
Based on the results above, the function space containing g0 of Eq. (8.10.18)
is included in W 1,qR (D; R) ⊂ X ′ which is smoother than the result of Theorem
8.5.2. From the fact that W 1,qR (D; R) ⊂ C 0 (D; R), it is thought that even
without applying the H 1 gradient method, a numerically unstable phenomenon
will not be generated. However, in order for the search vector ϑg to guarantee
C 0,1 class, the H 1 gradient method is required.

54

Chapter 8 Topology Optimization Problems of Density Variation Type8.10.4

Second-Order θ-Derivatives of Cost Functions

Furthermore, the second-order θ-derivatives of the cost functions of the mean
ﬂow resistance f0 and constraint function f1 with respect to the domain measureof ﬂow ﬁeld can be obtained. Based on the procedures looked at in Sect. 8.5.2,
let us also look here at obtaining the second-order θ derivativs of f0 and f1 .
Firstly, let us think about the second-order θ-derivative of f0 . With respect
to Hypothesis 8.5.4 (1), b is assumed not to be a function of θ. The relationshipcorresponding to Hypothesis 8.5.4 (2) is satisﬁed here.
The Lagrange function L0 of f0 is deﬁned in Eq. (8.10.11). Viewing (θ, u, p)
as a design variable, its admissible set and admissible direction is set as
S = { (θ, u, p) ∈ D × S × Q |LS (θ, u, p, v, q) = 0 for all (v, q) ∈ U × P } ,TS (θ, u, p) = { (ϑ, υ̂, π̂) ∈ X × U × P |LSθup (θ, u, p, v, q) [ϑ, υ̂, π̂] = 0 for all (v, π̂) ∈ U × P } .
The second-order θ-derivative of L0 with respect to arbitrary variations(ϑ1 , υ̂ 1 , π̂1 ) , (ϑ2 , υ̂ 2 , π̂2 ) ∈ TS (θ, u, p) becomes
L0(θ,u,p)(θ,u,p) (θ, u, p, v 0 , q0 ) [(ϑ1 , υ̂ 1 , π̂1 ) , (ϑ2 , υ̂ 2 , π̂2 )]= L0θθ (θ, u, p, v 0 , q0 ) [ϑ1 , ϑ2 ] + L0θup (θ, u, p, v 0 , q0 ) [ϑ1 , υ̂ 2 , π̂2 ]+ L0θup (θ, u, v 0 ) [ϑ2 , υ̂ 1 , π̂1 ] + L0upup (θ, u, v 0 ) [υ̂ 1 , π̂1 , υ̂ 2 , π̂2 ] .
(8.10.21)
Each term on the right-hand side of Eq. (8.10.21) becomes
L0θθ (θ, u, p, v 0 , q0 ) [ϑ1 , ϑ2 ]
Z n
o
2
=
ψ ′′ (ϕ (θ)) (ϕ′ (θ)) + ψ ′ (ϕ (θ)) ϕ′′ (θ) u · v 0 ϑ1 ϑ2 dx, (8.10.22)
DZ
L0θup (θ, u, p, v 0 , q0 ) [ϑ1 , υ̂ 2 , π̂2 ] =ψ ′ (ϕ (θ)) ϕ′ (θ) υ̂ 2 · v 0 ϑ1 dx,
DZ
L0θup (θ, u, p, v 0 , q0 ) [ϑ2 , υ̂ 1 , π̂1 ] =(8.10.23)
ψ ′ (ϕ (θ)) ϕ′ (θ) υ̂ 1 · v 0 ϑ2 dx,

D(8.10.24)
L0upup (θ, u, p, v 0 , q0 ) [υ̂ 1 , π̂1 , υ̂ 2 , π̂2 ] = 0.

(8.10.25)

Here, the fact that u − uD , v 0 − uD , υ̂ 1 and υ̂ 2 become 0Rd on ∂D was used.
In this case, with respect to arbitrary variation (ϑj , υ̂ j , π̂j ) ∈ TS (θ, u, p) forj ∈ {1, 2}, the Fréchet partial derivative of the Lagrange function LS of thestate determination problem establishes the equation
LSθup (θ, u, p, v, q) [ϑ, υ̂, π̂]
Z n



=
−µ ∇υ̂ ⊤ · ∇v ⊤ − ψ ′ (ϕ (θ)) ϕ′ (θ) u · v ϑD8.10 Topology Optimization Problem of Stokes Flow Field55o− ψ (ϕ (θ)) υ̂ · v + π̂∇ · v + q∇ · υ̂ dx=0(8.10.26)

with respect to an arbitrary (v, q) ∈ U × P . Here, the fact that v and υ̂ are
0Rd on ΓD as well as Eq. (8.10.5) were used.
Here, the next assumption is set up. At a local minimum point of a mean
ﬂow resistance minimization problem (Problem 8.10.3), it is thought that theﬂuid content ϕ (θ) converges to 0 and 1, and the terms introduced in order toprovide ﬂow resistance becomes suﬃciently small in actual ﬂow ﬁeld of ϕ (θ) ≈ 1.
Hence, in Eq. (8.10.26), it is assumed that
Z n


o

−µ ∇υ̂ ⊤· ∇v ⊤ + π̂j ∇ · v + q∇ · υ̂ j dx = 0(8.10.27)
jD

is established. In this case, the conditionυ̂ j = −ψ ′ (ϕ (θ)) ϕ′ (θ)ϑj uψ (ϕ (θ))

in D(8.10.28)

can be obtained. Hence, substituting (υ̂ j , π̂j ) into (υ̂ 1 , π̂1 ) in Eq. (8.10.24) and(υ̂ 2 , π̂2 ) in Eq. (8.10.23), then
L0θup (θ, u, v 0 ) [ϑ1 , υ̂ 2 , π̂2 ]= L0θup (θ, u, p, v 0 , q0 ) [ϑ2 , υ̂ 1 , π̂1 ]Z2(ψ ′ (ϕ (θ)) ϕ′ (θ))=−
u · v 0 ϑ1 ϑ2 dx
ψ (ϕ (θ))
D(8.10.29)

is obtained.
Summarizing the results above, by substituting Eq. (8.10.29) andEq. (8.10.22) into Eq. (8.10.21), the second-order θ-derivative of mean ﬂow
resistance f0 becomesZ 2(ψ ′ ϕ′ )2h0 (ϑ1 , ϑ2 ) =ψ ′′ (ϕ′ ) + ψ ′ ϕ′′ − 2
u · v 0 ϑ1 ϑ2 dx
ψDZ(8.10.30)=
ψ1 β (α, θ) u · v 0 ϑ1 ϑ2 dx.
D

In the above equation, ψ ′ (ϕ) becomes Eq. (8.10.19) andψ ′′ (ϕ) =2α (1 + α)3

(ϕ + α).

(8.10.31)

When Eq. (8.1.1) is used in ϕ (θ), ϕ′ (θ) and ϕ′′ (θ) are given by Eq. (8.5.7) andEq. (8.9.27), respectively. Moreover, if ϕ (θ) is given by Eq. (8.1.2), these are
given by Eq. (8.5.8) and Eq. (8.9.28), respectively. Figure 8.19 shows the graphof β (α, θ). From β (α, θ) < 0 and the self-adjoint relationship u · v 0 = u · u > 0,
it can be conﬁrmed that h0 ( · , · ) is not coercive.

56

Chapter 8 Topology Optimization Problems of Density Variation Type(a) ϕ (θ) = tan−1 θ/π + 1/2(b) ϕ (θ) = (tanh θ + 1) /2

Fig. 8.19: Coeﬃcient functions β (α, θ) at second-order θ-derivative of mean
ﬂow resistance.
On the other hand, the second-order θ-derivative of f1 (θ) becomes
Eq. (8.9.26). The graph of ϕ′′ (θ) is shown in Fig. 8.13.
In this way, the second-order θ-derivative of the cost function f0 in mean
ﬂow resistance minimization problem is not coercive and the second-orderθ-derivative of the constraint function f1 does not become coercive either.
Hence, if the Newton method (Problem 8.6.6) is to be used with respect to amean ﬂow resistance minimization problem, there is a need to use an appropriate
bilinear form aX (ϑgi , ψ) to capture coerciveness.

8.10.5

Second-Order θ-Derivative of Cost Function Using
Lagrange Multiplier Method

When the Lagrange multiplier method is used to obtain the second-orderθ-derivative of the mean ﬂow resistance f0 , it becomes as follows. Using thesame discussion as Section 7.5.4, we ﬁx ϑ1 and deﬁne the Lagrange function forf˜0′ (θ) [ϑ1 ] = ⟨g0 , ϑ1 ⟩ in Eq. (8.10.17) by
LI0 (θ, u, p, w0 , r0 ) = ⟨g0 , ϑ1 ⟩ − LS (θ, u, p, w0 , r0 ) ,(8.10.32)

where LS is given by Eq. (8.10.8), and (w0 , r0 ) ∈ U × P is the adjoint variable
provided for (u, p) in g 0 .
2
With respect to arbitrary variations (ϑ2 , û, p̂, ŵ0 , r̂0 ) ∈ X × (U × P ) of
(θ, u, p, w0 , r0 ), the Fréchet derivative of LI0 is written as
LI0′ (θ, u, p, w0 , r0 ) [ϑ2 , û, p̂, ŵ0 , r̂0 ]= LI0θ (θ, u, p, w0 , r0 ) [ϑ2 ] + LI0up (θ, u, p, w0 , r0 ) [û, p̂]+ LI0w0 r0 (θ, u, p, w0 , r0 ) [ŵ0 , r̂0 ] .

(8.10.33)

The third term on the right-hand side of Eq. (8.10.33) vanishes if (u, p) is thesolution of the state determination problem.
The second term on the right-hand side of Eq. (8.10.33) is
LI0up (θ, u, p, w0 , r0 ) [û, p̂]8.10 Topology Optimization Problem of Stokes Flow Field
Z n

 
⊤
=
2ψ ′ ϕ′ u · ûϑ1 + µ ∇w⊤+ ψ (ϕ (θ)) w0 · û
0 · ∇ûDo− p̂∇ · w0 − r0 ∇ · û dx.

57(8.10.34)

Here, the condition that Eq. (8.10.34) is zero for arbitrary (û, p̂) ∈ U × P isequivalent to setting (w0 , r0 ) to be the solution of the following adjoint problem.
Problem 8.10.4 (Adjoint problem of (w0 , r0 ) with respect to ⟨g0 , ϑ1 ⟩)Under the assumption of Problem 8.10.2, let ϑ1 ∈ X be given. Find
(w0 , r0 ) = (w0 (ϑ1 ) , r0 (ϑ1 )) ∈ U × P satisfying

⊤
⊤
′ ′ ⊤
−∇⊤ µ∇w⊤
0 + ψw 0 + ∇ r0 = −2ψ ϕ v 0 ϑ1
∇ · w0 = 0Zw 0 = 0Rd

in D,

in D,
on ∂D,

r0 dx = 0.
D□Finally, the ﬁrst term on the right-hand side of Eq. (8.10.33) becomes
LI0θ (θ, u, p, w0 , r0 ) [ϑ2 ]
Z n

o
2
=
ψ ′′ (ϕ′ ) + ψ ′ ϕ′′ u · uϑ1 + ψ ′ ϕ′ u · w0 (ϑ1 ) ϑ2 dx.
DHere, (u, p) and (w0 (ϑ1 ) , r0 (ϑ1 )) are assumed to be the weak solutions ofProblems 8.10.2 and 8.10.4, respectively. If we denote fi (θ, u, p) here by f˜i (θ),
we have the relation:
LI0θ (θ, u, p, w0 (ϑ1 ) , r0 ) [ϑ2 ] = f˜0′′ (θ) [ϑ1 , ϑ2 ]= ⟨gH0 (θ, ϑ1 ) , ϑ2 ⟩ ,(8.10.35)

where the Hesse gradient gH0 of the mean ﬂow resistance is given by2
gH0 (θ, ϑ1 ) = ψ ′′ (ϕ′ ) + ψ ′ ϕ′′ u · uϑ1 + ψ ′ ϕ′ u · w0 (ϑ1 ) .

(8.10.36)

If the same relation with Eq. (8.10.27) is satisﬁed in Problem 8.10.4,
w0 (ϑ1 ) = −ψ ′ ϕ′
v 0 ϑ1ψ(8.10.37)

holds. Substituting Eq. (8.10.37) into Eq. (8.10.36), it can be conﬁrmed thatEq. (8.10.35) accords with Eq. (8.10.30).

58

Chapter 8 Topology Optimization Problems of Density Variation Type(a) Initial density and boundary condition(b) H 1 gradient method (k = 100)(c) H 1 Newton method (k = 100)(d) H 1 Newton method(Hesse gradient, k = 100)

Fig. 8.20: Numerical example of mean ﬂow resistance minimization: densities(black and white are inverted).

(a) Initial density(b) H 1 Newton method

Fig. 8.21: Numerical example of mean ﬂow resistance minimization problem:
streamlines

Cost function

8.10 Topology Optimization Problem of Stokes Flow Field

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.0

0

20

40

60

80

100

(a) Cost functions

0.0

0

2

2

4

4

0

20

40

60

80

100

(c) Gradient of f0 on the search path

6

0

0.12

0.10

0.10

0.08

0.08

0.06

0.06

0.04

0.04

0.02

0.02
0

20

40

60

30

40

50

80

100

(e) Hessian of f0 on the search path10
20
Search distance

30

40

50

(d) Gradient of f0 on the search path(search distance)

0.12

0.00

10
20
Search distance(b) Cost functions (search distance)

0

6

0

59

0.00

0

10
20
Search distance

30

40

50

(f) Hessian of f0 on the search path(search distance)

Fig. 8.22: Numerical example of mean ﬂow resistance minimization: costfunctions and gradients and Hessians of f0 on the search path (gL : H 1 gradientmethod, hL , gL : H 1 Newton method, gH0 , h1 , gL : H 1 Newton method using
Hesse gradient).

60

Chapter 8 Topology Optimization Problems of Density Variation Type
3.0
2.5
3.0

2.0

2.5
1.5

2.0
1.5

1.0

2

1.0

1

0.5
20

0.5
40

60

80

100

(a) Iteration history

0.5

1

1.5

2

2.5

3

(b) (k − 1)-th vs. k-th plot

Fig. 8.23: Numerical example of mean ﬂow resistance minimization problem:
distance θ(k) − θ∗ X from an approximate minimum point θ∗ (gL : the gradientmethod, hL , gL : the Newton method，gH0 , h1 , gL : the Newton method using
the Hesse gradient).

8.10.6

Numerical Example

The results of mean ﬂow resistance minimization for a two-dimensional Stokes
ﬂow ﬁeld around an isolated object are shown in Figs. 8.20 to 8.23. The
boundary condition of a state determination problem is assumed on the outerboundary with a uniform ﬂow ﬁeld in the horizontal direction as shown in
Fig. 8.20 (a). With respect to the initial θ, a two-dimensional Gaussian
distribution was assumed in order for it to be an element of the admissible
set D. Also in this example, a domain in which the density is constrained was
not set. The programs were written using the programming language FreeFEM(https://freefem.org/) [16] using the ﬁnite element method. In the ﬁnite
element analyses of the Stokes problem, triangular elements of the second orderwith respect to the velocity and of the ﬁrst order with respect to the pressure
were used. In the case using the H 1 Newton method, the routine of the H 1
Newton method was started at kN = 20. For further details, we recommend thereaders to also examine the exact codes in the programs.2
Figure 8.20 (b) to (d) show the densities obtained by the three methods(H 1 gradient method using gL = g0 + λ1 g1 , H 1 Newton method using hL =h0 + λ1 h1 and gL , and the H 1 Newton method using gH0 , h1 , gL ). In Fig. 8.21,
the streamlines at the initial density and at the optimal density obtained fromthe H 1 Newton method are shown, respectively. The streamlines are deﬁned as
the contour lines of the ﬂow function ψ : Ω (ϕ) → R such that the ﬂow speed u⊤
is given by (∂ψ/∂x2 , −∂ψ/∂x1 ) .
The graphs in Fig. 8.22 illustrate the histories of the cost functions as wellas the gradients and Hessians of the objective function f0 on the search path
2 See Electronic supplementary material.

8.11 Summary

61

Pk−1
with respect to the iteration number k and the search distance i=0 ϑg(i) X
on X. In this ﬁgure, f0init denotes the value of f0 at the initial density.
Also, the value of c1 is taken to be the domain integral of the initial density(volume). The gradient of f0 on the search path was calculated using the
Lagrange function L = L0 + λ1 f1 as g0 + λ1 g1 , ϑg(k) / ϑg(k) X . On
the other hand, the Hessian of f0 on the search path was computed via2hL ϑg(k) , ϑg(k) / ϑg(k) X . In the case of the Newton method using the2
Hesse gradient, the ratio gH0 , ϑg(k) + λ1 h1 ϑg(k) , ψ / ϑg(k) X was used
to calculate the Hessian. The norm ϑg(i) X of the i-th search vector is deﬁned
by Eq. (8.9.35). The computational times until k = 100 by PC were 10.839,
14.763 and 17.046 sec by the H 1 gradient method, the H 1 Newton method andthe H 1 Newton method using the Hesse gradient, respectively.
Looking at the graphs in Fig. 8.22 (c), it also seems that the convergence
speed with respect to the iteration number k is faster when the H 1 Newton
method is applied than when the H 1 gradient method is used. However, thisincrease in convergence speed might have been due to the fact that cD1 and cD0in Eq. (8.6.7) are set with smaller values that had made the step sizes larger.
Meanwhile, we noticed that the search paths for the three methods are the same
as evident from the three coinciding graphs plotted in Fig. 8.22 (d). Moreover,from Fig. 8.22 (d) and (f), it can be observed that the point of convergence isa local minimum and that both the ﬁrst derivative and the Hessian eventuallydiminish to zero. The reason behind these observed behaviors of the methods
is the same as the ones given at the end of Sect. 8.9.6. Moreover, in the case of
the Stokes ﬂow ﬁeld, although h0 ( · , · ) itself is not coercive as shown in Sect.
8.10.4, it can be conﬁrmed that the point of convergence is a minimum point,since the Hessian of f0 on the search path is positive valued.
In addition, Fig. 8.23 (a) shows the graphs of the distance θ(k) − θ∗ X froman approximate minimum point θ∗ obtained by the three methods with respect
to the iteration number k. The approximate minimum point θ∗ was substituted
with the numerical solution of θ when the iteration time is taken larger than thegiven value in the H 1 Newton method. From this ﬁgure, it can be conﬁrmedthat the convergence orders for the results via the H 1 Newton methods arehigher compared to that of the ﬁrst order. However, from Fig. 8.23 (b), plotting
the k-th distance θ(k) − θ∗ X with respect to the (k − 1)-th distance, theconvergence order of the H 1 Newton method is less than the second order but
is greater than the ﬁrst order. The reason behind this result is considered thesame as that stated at the end of Sect. 8.9.6.

8.11

Summary

In Chap. 8, the problem for seeking optimal hole positions with respect toa domain, in which a boundary value problem of partial diﬀerential equation
is deﬁned, is constructed as topology optimization problem of θ-type and its
solution was looked at in detail. The key points are as below:

62

Chapter 8 Topology Optimization Problems of Density Variation Type(1) If the characteristic function of a domain is chosen as a design variable,
it is known that regularity is insuﬃcient and an optimization problem
cannot be constructed (beginning of Chap. 8).
(2) If a density is chosen as the design variable, a topology optimization
problem can be constructed (beginning of Chap. 8). However, the setof functions whose range is limited to [0, 1] does not become a linear
space. Hence, by choosing a function θ ∈ X = H 1 (D; R) whose range
is not constrained as the design variable and providing the density asa sigmoid function of θ, a topology optimization problem of θ-type canbe constructed based on the framework of the abstract optimal design
problem shown in Chap. 7 (Sect. 8.1).
(3) When a Poisson problem is chosen as a state determination problem (Sect.
8.2), a topology optimization problem of θ-type will be constructed as
Problem 8.3.2 (Sect. 8.3).
(4) The θ-derivatives of cost functions can be obtained by the Lagrange
multiplier method (Sect. 8.5.1). However, such θ-derivatives are notnecessarily going to be in X (Remark 8.5.3). Moreover, the second-orderθ-derivatives of cost functions can be sought by substituting in the
θ-derivative of the solution of a state determination problem into thesecond-order θ-derivative of the Lagrange function (Sect. 8.5.2).
(5) The descent vectors of cost functions can be obtained using the
θ-derivatives of cost functions by a gradient method (H 1 gradient method)
on X = H 1 (D; R) (Sect. 8.6.1). The solution of the H 1 gradient method
is included in the admissible set apart from the singular points (Theorem
8.6.5). Furthermore, if the second-order θ-derivatives of cost functions canbe calculated, the descent vectors of the cost functions can be sought via
the H 1 Newton method (Sect. 8.6.2).
(6) The solution to the topology optimization problem of θ-type can beconstructed using the same framework as the gradient method with respect
to a constrained problem and the Newton method with respect to a
constrained problem shown in Chap. 3 (Sect. 8.7.1 and Sect. 8.7.2).
(7) When the numerical solutions of a state determination problem, adjoint
problem, and H 1 gradient method are to be sought via the ﬁnite element
method, and a ﬁrst-order ﬁnite element is used to seek the search vector
ϑg , the error of the ﬁnite element solution reduces linearly with respect tothe maximum diameter of the ﬁnite element (Theorem 8.8.5).
(8) When a linear elastic problem is taken to be a state determination
problem, and a mean compliance and a function for domain measureare chosen as object and constraint cost functions, the θ-derivatives andsecond-order θ-derivatives of the cost functions can be obtained (Sect.
8.9.3).

8.12 Practice Problems

63(9) If the Stokes problem is taken to be a state determination problem, anda mean ﬂow resistance minimization problem under a domain measureconstraint is considered, the θ-derivatives and second-order θ-derivativesof the cost functions can be obtained (Sect. 8.10.3).

8.12

Practice Problems

8.1 When a θ-type Poisson problem (Problem 8.2.3) is made into a state
determination problem, what would be the cost function such that the
self-adjoint relationship holds? Moreover, show its θ-derivative.
8.2 Change the extended Poisson problem deﬁned in Chap. 5 (Problem
5.1.3) to θ-type, and formulate a topology optimization problem of θ-typeusing the extended Poisson problem as a state determination problem, theobject cost function such that the self-adjoint relationship holds, and theconstraint cost function with respect to the domain measure. Moreover,show the KKT conditions with respect to that problem.
8.3 When many cost functions are deﬁned and the maximum values of them
have to be minimized, the β method is known as a way to construct an
optimal design problem [39]. If the topology optimization problem ofθ-type is rewritten with the β method, we obtain the following:
Problem 8.12.1 (The β method) Let D and S be Eq. (8.1.4) andEq. (8.2.2), respectively, and f1 , . . . , fm : X × S → R be given as
Eq. (8.3.1). Moreover, let β ∈ R. In this case, obtain θ which satisﬁes:min(θ,u−uD )∈D×S{ β | f1 (θ, u) ≤ β, . . . , fm (θ, u) ≤ β, Problem 8.2.3} .
□Show the KKT conditions with respect to this problem. Moreover, show
the method to determine the Lagrange multipliers when solving thisproblem using the H 1 gradient method (Sect. 8.7.1) with respect to a
constrained problem.
(Supplement) The reason that the β method is preferred is shown as
follows. Even if there are many cost functions, the Lagrange multiplierswith respect to cost functions for which inequality constraints are notactive are zero. Hence, there are no reasons to seek their θ-derivatives.
8.4 In Practice 1.2, the gradient g 0 of mean compliance f0 with respect
to variation of cross-sectional area was sought by using the gradient of
maximization problem of the potential energy with respect to variationof a and the minimum condition of the potential energy with respect
to variation of u. With respect to the mean compliance minimization

64

Chapter 8 Topology Optimization Problems of Density Variation Typeproblem (Problem 8.9.3) of a d ∈ {2, 3}-dimensional linear elastic body,think about the problem using

Z 
1 α
π (θ, u) =ϕ (θ) S (u) · E (u) − b (θ) · u dx
2
D
Z
Z
−
pN · u dγ −(u − uD ) · (ϕα (θ) S (v 0 ) ν) dγΓN

ΓDas the potential energy and seeking (θ, u) that satisﬁes
max min π (θ, u) .
θ∈D u∈U

In this case, show that the θ gradient when using u satisfying minu∈U π
is the same as half of g0 in Eq. (8.9.14).
8.5 In Practice 1.8, the gradient g 0 of mean ﬂow resistance f0 with respect
to variation of cross-sectional area was obtained using the gradient ofminimization problem of a formal potential energy for a dissipative
system with respect to variation of a and maximum condition of theformal potential energy with respect to variation of p. With respect
to the mean ﬂow resistance minimization problem (Problem 8.10.3) ofd ∈ {2, 3}-dimensional Stokes ﬂow ﬁeld,
Z n

 1
1
π (θ, u, p) =µ ∇u⊤ · ∇u⊤ + ψ (ϕ (θ)) u · u − p∇ · u
2
D 2
Z
o
− b · u dx −(u − uD ) · (µ∂ν u − pν) dγ∂Dis used to seek (θ, u, p) which satisﬁes
min min max π (θ, u, p) .
θ∈D u∈U p∈P

In this case, show that the θ gradient of π when using (u, p) satisfying
minu∈U maxp∈P π is the same as half of g0 in Eq. (8.10.17).

References[1] Aage, N., Poulsen, T. H., Gersborg-Hansen, A., and Sigmund, O.
Topology optimization of large scale Stokes ﬂow problems. Structural andMultidisciplinary Optimization, 35:175–180, 2008.
[2] Allaire, G. Shape Optimization by the Homogenization Method. Springer,
New York, 2002.
[3] Avellaneda, M. Optimal bounds and microgeometries for elastic two-phasecomposites. SIAM J. Appl. Math., 47:1216–1228, 1987.
[4] Azegami, H., Kaizu, S., and Takeuchi, K. Regular solution to topology
optimization problems of continua. JSIAM Letters, 3:1–4, 2011.
[5] Bendsøe, M. P. Optimal shape design as a material distribution problem.
Structural Optimization, 1:193–202, 1989.
[6] Bendsøe, M. P. and Kikuchi, N. Generating optimal topologies in structuraldesign using a homogenization method. Computer Methods in Applied
Mechanics and Engineering, 71:197–224, 1988.
[7] Bendsøe, M. P. and Sigmund, O. Topology Optimization: Theory, Methods
and Applications. Springer, Berlin; Tokyo, 2003.
[8] Borrvall, T. and Petersson, J. Topology optimization of ﬂuids in Stokes
ﬂow. International Journal for Numerical Methods in Fluids, 41:77–107,
2003.
[9] Diaz, A. R. and Kikuchi, N. Solutions to shape and topology eigenvalue
optimization problems using a homogenization method. International
Journal for Numerical Methods in Engineering, 35:1487–1502, 1992.
[10] Diaz, A. R. and Sigmund, O. Checkerboard patterns in layout optimization.
Structural Optimization, 10:40–45, 1995.
[11] Evgrafov, A. Topology optimization of slightly compressible ﬂuids. Journalof Applied Mathematics and Mechanics / Zeitschrift f´’ur Angewandte
Mathematik und Mechanik, 86:46–62, 2006.
65

66

Chapter 8 Topology Optimization Problems of Density Variation Type[12] Gersborg-Hansen, A., Sigmund, O., and Haber, R. B.
Topology
optimization of channel ﬂow problems. Structural and Multidisciplinary
Optimization, 3:181–192, 2005.
[13] Guest, J. K. and Prévost, J. H. Topology optimization of creeping
ﬂuid ﬂows using a Darcy-Stokes ﬁnite element. International Journal for
Numerical Methods in Engineering, 66:461–484, 2006.
[14] Guest, J. K. and Prévost, J. H. Design of maximum permeability material
structures. Computer Methods in Applied Mechanics and Engineering,
196:1006–1017, 2007.
[15] Haslinger, J. and Mäkinen, R. A. E. Introduction to Shape Optimization:
Theory, Approximation, and Computation. SIAM, Philadelphia, 2003.
[16] Hecht, F. New development in freefem++.
Mathematics, 20(3-4):251–265, 2012.

Journal of Numerical[17] Kawamoto, A., Matsumori, T., Yamasaki, S., Nomura, T., Kondoh,T., and Nishiwaki, S. Heaviside projection based topology optimization
by a PDE-ﬁltered scalar function.
Structural and Multidisciplinary
Optimization, 44:19–24, 2011.
[18] Kohn, R. V. and Lipton, R. Optimal bounds for the eﬀective energy of a
mixture of isotropic, incompressible, elastic materials. Archive for Rational
Mechanics and Analysis, 102(4):331–350, 1988.
[19] Kohn, R. V. and Strang, G. Optimal design and relaxation of variational
problems, part 1. Communications on Pure and Applied Mathematics,
39:1–25, 1986.
[20] Kohn, R. V. and Strang, G. Optimal design and relaxation of variational
problems, part 2. Communications on Pure and Applied Mathematics,
39:139–182, 1986.
[21] Kohn, R. V. and Strang, G. Optimal design and relaxation of variational
problems, part 3. Communications on Pure and Applied Mathematics,
39:353–377, 1986.
[22] Kohn, R. V. and Strang, G. Optimal design in elasticity and plasticity.
International Journal for Numerical Methods in Engineering, 22:183–188,
1986.
[23] Lazarov, B. S. and Sigmund, O. Filters in topology optimization based on
Helmholtz-type diﬀerential equations. International Journal for Numerical
Methods in Engineering, 86(6):765–781, 2011.
[24] Lurie, K. A., Cherkaev, A. V., and Fedorov, A. V. Regularization of optimal
design problems for bars and plates, part 1. Journal of Optimization Theory
and Applications, 37:499–522, 1982.

8.12 Practice Problems

67[25] Lurie, K. A., Cherkaev, A. V., and Fedorov, A. V. Regularization of optimal
design problems for bars and plates, part 2. Journal of Optimization Theory
and Applications, 37:523–543, 1982.
[26] Ma, Z. D., Kikuchi, N., and Hagiwara, I. Structural topology and shape
optimization for a frequency response problem. Computational Mechanics,
13:157–174, 1993.
[27] Matsui, K. and Terada, K. Continuous approximation of material
distribution for topology optimization. International Journal for Numerical
Methods in Engineering, 59:1925–1944, 2004.
[28] Mlejnek, H. P. and Schirrmacher, R. An engineer’s approach to optimal
material distribution and shape ﬁnding. Computer Methods in Applied
Mechanics and Engineering, 106:1–26, 1993.
[29] Murai, D. Error analyses in numerical solutions for initial and boundary
value problems of partial diﬀerential equations and shape optimization
problems (in Japanese). PhD thesis, Nagoya University, 2011.
[30] Murai, D. and Azegami, H. Error analysis of H1 gradient method for
topology optimization problems of continua. JSIAM Letters, 3:73–76, 2011.
[31] Murat, F. Contre-exemples pour divers problémes ou le contrôle intervient
dans les coeﬃcients. Annali di Matematica Pura ed Applicata, Serie 4,
112:49–68, 1977.
[32] Nishiwaki, S., Izui, K., and Kikuchi, N.
Japanese). Maruzen, Tokyo, 2013.

Topology Optimization (in[33] Olesen, L. H., Okkels, F., and Bruus, H.
A high-level
programming-language implementation of topology optimization applied
to steady-state Navier-Stokes ﬂow. International Journal for Numerical
Methods in Engineering, 65:975–1001, 2006.
[34] Rahmatalla, S. F. and Swan, C. C. A Q4/Q4 continuum structural
topology optimization implementation. Structural and Multidisciplinary
Optimization, 27:130–135, 2004.
[35] Rozvany, G., Zhou, M., and Birker, T. Generalized shape optimization
without homogenization. Structural Optimization, 4:250–254, 1992.
[36] Sigmund, O. and Maute, K. Topology optimization approaches: A
comparative review.
Structural and Multidisciplinary Optimization,
48:1031–1055, 2013.
[37] Sigmund, O. and Petersson, J. Numerical instabilities in topology
optimization: A survey on procedures dealing with checkerboards,
mesh-dependencies and local minima. Structural Optimization, 16:68–75,
1998.

68

Chapter 8 Topology Optimization Problems of Density Variation Type[38] Suzuki, K. and Kikuchi, N. A homogenization method for shape andtopology optimization. Computer Methods in Applied Mechanics andEngineering, 93:291–318, 1991.
[39] Taylor, J. E. and Bendsøe, M. P. An interpretation for min-max structuraldesign problems including a method for relaxing constraints. International
Journal of Solids and Structures, 20:301–314, 1984.
[40] Wang, F., Lazarov, B. S., and Sigmund, O. On projection methods,
convergence and robust formulations in topology optimization. Structural
and Multidisciplinary Optimization, 43:767–784, 2011.
[41] Yang, R. J. and Chuang, C. H. Optimal topology design using linear
programming. Computers & Structures, 52:265–275, 1994.