Research Reports of National Institute of Technology, Nagaoka College, Vol. 59, 2024

Paper

Comparing the topology optimization resultsof the stiffest structureusing two types of topological derivatives
Takahiko KURAHASHI1, Manaya NAKAZAWA2, Masayuki KISHIDA3,4,
Kousuke SHIMIZU5 and Tetsuro IYAMA6

1

Department of Mechanical Engineering, Institute of GIGAKU, Nagaoka University of Technology
Mechanical Engineering, Bachelor’s Program, School of Engineering, Nagaoka University of Technology

2(Current affiliation : ADTEC Engineering Co.,Ltd.)
3

Department of Mechanical Engineering, National Institute of Technology, Gifu College

4

Science of Technology Innovation, 5-year Integrated Doctoral Program, Graduate School of Engineering,
Nagaoka University of Technology

Mechanical Engineering, Master’s Program, Graduate School of Engineering, Nagaoka University of Technology

56Department of Mechanical Engineering, National Institute of Technology, Nagaoka College

This study examined the effect of varying the approach to obtain the topological derivative and presents thenumerical results of level set-based topology optimization for a maximally stiff structure problem. To perform a
topology optimization analysis, the performance function was first defined by the strain energy of the structure.
The problem was to determine the optimal topology to minimize the performance function under constraint
conditions, that is, the governing equation and boundary condition. The adjoint variable method was introduced toaddress the minimization problem of the performance function under constraint conditions. The optimal topologyof the structure was obtained by updating the level-set function, which was achieved by solving the reactiondiffusion equation. The reaction term of the reaction-diffusion equation was expressed by the topologicalderivative, that is, the gradient of the performance function extended by the adjoint variable and the governing
equation with respect to the level-set function. In this study, we varied the method to obtain the topologicalderivative in level-set-based topology optimization and performed numerical experiments. The finite element
method was applied to solve the structural deformation problem.
Key Words: Level-set-based topology optimization, Maximally stiff structure, Adjoint variable method,
Topological derivative, Finite element method

1. Introduction

resemble the tissue structure of the bones of livingorganisms1). The present study focuses on level-setbased topology optimization, that is, a solution
methodalongside
density-based
topology
optimization. In some topological optimization
problems, the induced topological derivative renders

Since three-dimensional (3D) printers have gained
prevalence, topology optimization analysis has
become a focus of attention. Recently, it has beenreported that the results of topology optimization

20

Comparing the topology optimization results of the stiffest structure using two types of topological derivatives
it difficult to compute the optimal topology2). The
topological derivative is specified by the gradient of
the objective function with respect to the level-set
function, but the convergence rate of the objectivefunction and the accuracy of the optimal solution
depends on the gradient, that is, the search direction,such as the relationship between gradient- andNewton-based methods. Therefore, numericalexperiments of topology optimization are conductedin this study to maximize stiffness by varying the
topological derivative.

symbol indicates the specified value.
𝑲𝒖 = 𝒇̂ on Γ1
𝒖=𝒖{
𝒇 = 𝒇̂ on Γ2

1

1

(3)

Here, 𝜎𝑖𝑗 , 𝜖𝑖𝑗 , 𝐷𝑖𝑗𝑘𝑙 and 𝑢𝑖 indicate components ofstress and strain tensors, the elasticity coefficienttensor and the displacement components,
respectively. If the finite element Galerkin procedure
is applied to discretize Eqs. (1)–(3), then the finite
element equation shown in Eq. (4) can be derived3).
𝒖𝑒 and 𝒇𝑒 denote displacement and external force
vectors, respectively. 𝑲𝑒 represents the stiffness
matrix and is calculated using Eq. (5). 𝑩𝑒 and 𝑫𝑒
indicate the matrix B and the elasticity coefficientmatrix, respectively.
(4)

𝑲𝑒 = ∫Ω 𝑩𝑇𝑒 𝑫𝑒 𝑩𝑒 𝑑Ω

(5)

𝑒

(8)

The problem is to determine the optimal topology tominimize the performance function 𝐽. To obtain thedisplacement vector, the finite element equation
shown in Eq. (6) must be solved. Therefore, this is aminimization problem with a constraint condition, as
shown in Eq. (6). Using the adjoint variable method,
the performance function is extended by the adjoint
variable vector 𝝀 and finite element equation, Eq. (6),and the extended performance function is expressedas shown in Eq. (9). Eq. (9) is referred to as the
Lagrange function. Here, the stiffness matrix 𝑲 isexpressed by the characteristic function 𝜒, and thecharacteristic function 𝜒 is the function of the levelset function 𝜙. If the level-set function 𝜙 is positive,then the element is in the material domain. If thelevel-set function 𝜙 is negative, then the element isin the void domain. If the level-set function 𝜙 is zero,then the boundary is between the material and void
domains. The characteristic function 𝜒 is 1 when theelement is in the material domain, and 0 when theelement is in the void domain. The stiffness matrix 𝑲
is employed when 𝜒 is 1. The stiffness matrix 𝑲 isexpressed by a zero matrix when 𝜒 is 0.

(1)
(2)

𝑲𝑒 𝒖𝑒 = 𝒇𝑒

1

𝐽 = 2 𝒖𝑇 𝒇 = 2 𝒖𝑇 𝑲𝒖

The equivalent, stress–strain relation, and strain–
displacement relation equations are expressed as Eqs.
(1)–(3), respectively. These equations are written
using the summation convention.

𝜖𝑖𝑗 = 2 (𝑢𝑖,𝑗 + 𝑢𝑗,𝑖 )(7)Here, we consider the maximally stiff structuralproblem. The maximally stiff structure is a structurethat minimizes strain energy. Therefore, theperformance function is defined by the strain energy,
as shown in Eq. (8). Substituting Eq. (6) into theexternal force vector 𝒇 yields Eq. (8), which isrepresented by the stiffness matrix 𝑲 anddisplacement 𝒖.

2. Derivation of stationary condition ofextended performance function

𝜎𝑖𝑗,𝑗 = 0
𝜎𝑖𝑗 = 𝐷𝑖𝑗𝑘𝑙 𝜖𝑘𝑙(6)Superposing Eq. (4) for all elements, the finite
element equation for the entire domain is obtained asshown in Eq. (6). The boundary condition is definedby Eq. (7). Γ1 and Γ2 indicate the Dirichlet andNeumann boundaries, respectively, and the hat

1𝐽∗ = (𝒖𝑇 𝑲(𝜒(𝜙))𝒖 + 𝝀𝑇 (𝑲(𝜒(𝜙))𝒖 − 𝒇))2(9)

To control the complexity of the structure, a
regularization term is generally included in the21Takahiko KURAHASHI, Manaya NAKAZAWA, Masayuki KISHIDA, Kousuke SHIMIZU
and Tetsuro IYAMA
Lagrange function in level-set-based topology
relationship.
optimization. To derive the discretized regularization
term, we introduce the steady-state diffusion
𝜕𝐽∗∗= 𝑲(𝜒(𝜙))𝒖 − 𝒇 = 𝟎
𝜕𝝀equation shown in Eq. (10). 𝜏 and 𝜙 represent theregularization parameter and level-set function,
𝜕𝐽∗∗= 𝑲(𝜒(𝜙))𝒖 + 𝝀𝑇 𝑲(𝜒(𝜙))respectively. The finite element Galerkin procedure𝜕𝒖
is applied to Eq. (10), which results in Eq. (11). 𝑯= 𝒇 + 𝑲𝑇 (𝜒(𝜙))𝝀 = 𝟎
and 𝒒 indicate the diffusion matrix and heat-flux
𝑲(𝜒(𝜙))𝒖 = 𝒇vector, respectively. The boundary conditions for Eq.
𝑲𝑇 (𝜒(𝜙))𝝀 = −𝒇(11) are shown in Eq. (12).
𝒖 = −𝝀𝜏𝜙,𝑖𝑖 = 0
𝜏𝑯𝝓 = 𝒒̂ on Γ1
𝝓=𝝓{
̂ on Γ2
𝒒=𝒒

(14)

(15)
(16)
(17)
(18)

In addition, the level-set function vector 𝝓 is updated
by the gradient of the modified Lagrange function 𝐽∗∗
with respect to the level-set function vector 𝝓. The
update of the level-set function is performed using
the Allen–Cahn equation, which is discretized by thefinite element method (see Eq. (19)). Here, 𝑴𝑒

(10)
(11)
(12)

The regularization term, expressed in Eq. (11) is
added to the Lagrange function, 𝐽∗ . In addition, a
volume constraint condition is added to the Lagrangefunction 𝐽∗ . Finally, the modified Lagrange function
𝐾 is defined as shown in Eq. (13). Λ, 𝑣𝑒 , and 𝑉𝑚𝑎𝑥
indicate the adjoint variable, element volume, andtarget volume, respectively. 𝜒𝑒 represents thecharacteristic function. The value of 𝜒𝑒 is 1when theelement is in the material domain. However the valueof 𝜒𝑒 is 0 when the element is in a defect domain, that
is, a void domain. 𝑚𝑥 refers to the number of
elements.

𝜕𝐽∗∗

indicates the mass matrix. If 𝜕𝝓 is calculated, thenEq. (19) can be represented by Eq. (20). Eq. (20) canbe expressed as Eq. (21), which shows the reactiondiffusion equation. Parameter 𝐶 is the normalization𝑚𝑥

𝜕𝐽∗ . Thisparameter and is expressed as 𝐶 = ∑𝑚𝑥
|
𝑒=1|
𝜕𝝓𝑒

parameter is adjusted for each iteration step.
𝜕𝐽∗∗

𝜕𝝓

𝑴𝑒 𝜕𝑡̃𝑒 = −𝛼 𝜕𝝓

(19)

𝑒

𝜕𝐽∗

𝜕𝝓

𝑴𝑒 𝜕𝑡̃𝑒 = −𝛼 (𝐶 𝜕𝝓 + τ𝑯𝑒 𝝓𝑒 + Λ)𝑒1𝐽∗∗ = 𝐽∗ + 𝜏𝝓𝑇 𝑯𝝓 + Λ(∑𝑚𝑥
𝑒=1 𝜒𝑒 (𝜙)𝑣𝑒 − 𝑉𝑚𝑎𝑥 ) (13)
2

𝜕𝐽∗

𝜕𝝓

𝑴𝑒 𝜕𝑡̃𝑒 + ατ𝑯𝑒 𝝓𝑒 = − (𝐶 𝜕𝝓 + Λ)𝑒

To obtain the stationary condition of the modified
Lagrange function 𝐽∗∗ , the first variation of themodified Lagrange function 𝐽∗∗ is calculated. Fromthe gradient of the modified Lagrange function 𝐽∗∗
with respect to the adjoint variable vector 𝝀, the finite
element equation for elastic deformation is obtained(see Eq. (14).). From the gradient of the modified
Lagrange function 𝐽∗∗ with respect to thedisplacement vector 𝒖 , the adjoint equation is
obtained (see Eq.(15).). Eq. (16) can be obtained
from Eq. (14), and Eq. (17) can be obtained from Eq.
(15). Comparing Eqs. (16) and (17), the relation
equation between the displacement vector and adjoint
variable vector is derived as shown in Eq. (18). This
relationship is referred to as the self-adjoint

(20)
(21)

𝜕𝐽∗

According to a previous study4), 𝜕𝝓𝑒 is calculatedusing Eq. (22), and 𝑲′𝑒 is represented by Eq. (23). In
matrix 𝑲′𝑒 , matrix 𝑨𝑒 in the two-dimensional planestress state is expressed as Eq. (24). Here, 𝐴1 and 𝐴2
are calculated using Eq. (25).
𝜕𝐽∗= 𝝀𝑇𝑒 𝑲′𝑒 𝒖𝑒 = −𝒖𝑇𝑒 𝑲′𝑒 𝒖𝑒
𝜕𝝓𝑒

(22)

𝑲′𝑒 = ∫Ω 𝑩𝑇𝑒 𝑨𝑒 𝑩𝑒 𝑑Ω

(23)

𝑒

𝑨𝑒 = [

22

𝐴1 + 2𝐴2
𝐴1
0𝐴1
𝐴1 + 2𝐴2
0

0
0]
𝐴2

(24)

Comparing the topology optimization results of the stiffest structure using two types of topological derivatives
matrix 𝑫𝑒 was employed instead of matrix 𝑨𝑒 in thecalculation of the topological derivative6). Therefore,
in the setting of matrix 𝑲′𝑒 , the matrix shown in Eq.
(23) was employed in case 1, the elasticity coefficientmatrix 𝑫𝑒 was applied to calculate the matrix 𝑲′𝑒
instead of matrix 𝑨𝑒 in case 2, and the numericalresults were compared.
The numerical results are as follows. Figs. 3 and 4
show the variations in the level-set and characteristic
functions, respectively. The optimal topology isexpressed using a characteristic function. If the levelset function is positive, then the characteristic
function is 1 and the element represents the material
domain. However, if the level-set function is negative,then the characteristic function is zero and theelement represents the void domain. Based on acomparison between cases 1 and 2, although aslightly different shape is obtained in the midst of theiterative computation, a similar shape is obtained in
the final iteration.

3(1−𝜈)(1−14𝜈+15𝜈 2 ){

𝐴1 = − 2(1+𝜈)(7−5𝜈)(1−2𝜈2 ) 𝐸
15(1−𝜈))(25)

𝐴2 = − 2(1+𝜈)(7−5𝜈) 𝐸𝜕𝐽∗
𝜕𝝓𝑒 is the topological derivative and is derived foreach topological-optimization problem. Matrix 𝑲′𝑒 ,that is, Eq. (23), is similar to stiffness matrix 𝑲𝑒
shown in Eq. (5). Therefore, in this study, numericalexperiments were performed by changing matrix
𝜕𝐽∗𝑲′𝑒 to stiffness matrix 𝑲𝑒 when computing 𝜕𝝓𝑒 . An
image diagram of the flow of the level-set-based
topology optimization is shown in Fig. 1 (seereference5) for details.).
Next
Nextiteration
iteration

Void
Void

Y
Material
Material

X

1.0 m

0−2

Fig.1 Image diagram of the flow of level-set-based topology

0

optimization.

0

07

2.0

3. Numerical experiments

Fig.2 Computational model

Numerical experiments for the level-set-based
topology optimization analysis of the maximally stiff
structure problem were conducted by varying themethod to obtain the topological derivative. The
computational model is shown in Fig. 2, and thecomputational conditions are listed in Table 1.
Deformation analysis was performed under a twodimensional plane-stress state.
The computational model is illustrated in Fig. 2. The
aim of this study was to determine the optimal
topology to minimize the strain energy of thestructure. In this study, the numerical results werecompared by changing the method for providing the
topological derivative. The elasticity coefficient

Table 1 Computational conditions
Number of nodes / elements

23

Mesh size[ ]Young’s modulus[

12800 / 13041
1.25×10-2]

206

Poisson’s ratio

0.3

Virtual time

0.7

Volume reduction rate [ ]55Number of steps for volume
reduction

100

Design domain for topology
optimization [ ]2.0×1.0C

1.0×10-6

Regularization parameter3.0×10-4

Takahiko KURAHASHI, Manaya NAKAZAWA, Masayuki KISHIDA, Kousuke SHIMIZU
and Tetsuro IYAMA
Number of
iteration

1

50

70

80

100

120

140

160

Case1 :A Matrix

Case2 :D Matrix

Number of
iteration

Case1 :A Matrix

Case2 :D Matrix

Fig. 3 Variation in level-set function.

Number of
iteration

1

50

70

80

1
0

Case1 :A Matrix

Case2 :D Matrix

Number of
iteration

100

120

140

160

1
0

Case1 :A Matrix

Case2 :D Matrix

Fig. 4 Variation in characteristic function.

24

Comparing the topology optimization results of the stiffest structure using two types of topological derivatives

4. Conclusions

topology optimization method based on the level set method
incorporating a fictitious interface energy, Computer

In this study, numerical results were obtained bychanging the method for obtaining the topologicalderivative in level-set-based topology optimization.
The maximally stiff structural problem was solved
for the optimization problem, and the cantilever
beam model was introduced as a numerical example.
The results showed that even if matrix 𝑨𝑒 in the
topological derivative can be replaced by anelasticity coefficient matrix 𝑫𝑒 , an appropriate
optimized topology can be obtained.
Additionally, we developed a new updated
equation for the design variable7). In the future, we
aim to apply the new update equation to the present
method and investigate the effect of theconvergence rate of the performance function.

Methods in Applied Mechanics and Engineering, 199,
pp.2876–2891, 2010.
7)

optimization analysis using a new acceleration gradientmethod based on the Taylor expansion and conjugate
direction, Journal of Fluid Science and Technology, Vol.17,
No.4, pp.1-12, 2022.
(Received July 19, 2023)

Acknowledgments: The computations were mainly
performed using computer facilities at Kyushu
University’s Research Institute for Information
Technology. We wish to thank the staff at Kyushu
University’s Research Institute for Information
Technology. We would like to thank Editage(www.editage.com) for English language editing.
References
1)

N. Aage, E. Andreassen, B.S. Lazarov and O. Sigmund:Giga-voxel computational morphogenesis for structuraldesign, Nature, Vol. 550, pp.84-86, 2017.

2)

T. Nakazono and H. Yoshikawa: An identification of defects
using topological derivatives for nondestructive evaluationof elastic materials with laser interferometry data,
Transactions of the Japan Society for Computational
Methods in Engineering, Vol. 22. No.1, pp.171-177, 2022.
(In Japanese.)3)

M. Kawahara : Finite element methods in incompressible,
adiabatic, and compressible flows, Springer, 2016.

4)

M. Otomori, T. Yamada, K. Izui and S. Nishiwaki: Matlab
code for a level set-based topology optimization method
using a reaction diffusion equation, Structural andMultidisciplinary Optimization, Vol.51, No.5, pp.11591172, 2014.

5)

S. Nishiwaki, K. Izui and N. Kikuchi: Topology
optimization, Maruzen Publishing Co.,Ltd., 2013. (InJapanese.)6)

H. Arata, M. Kishida and T. Kurahashi : Texture shape

T. Yamada, K. Izui, S. Nishiwaki and A. Takezawa : A

25