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

Paper

Density-based topology optimization
for maximally stiff structure problem
using Two-phase materials
Takahiko KURAHASHI1, Sakura ICHIKAWA2, Masayuki KISHIDA3,4,
Mizuki IKARASHI5 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 : OILES CORPORATION)
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

Herein, we present numerical results obtained after performing density-based topology optimization fora maximally stiff structure problem using two-phase materials. A maximally stiff structure was modeled tominimize the strain energy of the structure. The Young’s modulus expressed based on a two-phase material,
such as a functionally graded material at the junction point, was used in the density-based topology
optimization. The two-phase material was completely separated in the optimized structure, and the strainenergies before and after the separation were compared.
Key Words: Density-based topology optimization, Two-phase material, Maximally stiff structure, Adjoint
variable method, Finite element method

1. Introduction

manufacturing, separating the material used from thejunction point is desirable.
In this study, a topology optimization analysis for a
maximally stiff structure is performed, and the valueof the performance function is compared before andafter the separation of two-phase materials. The
performance function is defined by the strain energy,
and numerical experiments are performed using theMesserschmitt–Bölkow–Blohm (MBB) model.

In recent years, researchers have attempted to use
various materials to reduce the weight of automobile
components1). Additionally, researchers have begun
to use topology optimization technology to
manufacture automobile components2). Considering
these trends, the topology optimization of two-phasematerials3) is investigated in this study. In topology
optimization analysis using a two-phase material, theYoung’s modulus is defined by a two-phase material.
Owing to topology optimization, a structure in which
the Young’s modulus varies gradually at the junction
point is obtained. However, in automobile

2. Formulations for density-based topology
optimization problem

1

Takahiko KURAHASHI, Sakura ICHIKAWA, Masayuki KISHIDA, Mizuki IKARASHI
and Tetsuro IYAMA
The governing equation for the linear elastic
known as the Lagrange function.
deformation is expressed in Eqs. (1)–(3). Theseequations are written using summation convention.
1∗𝑇
𝑇

𝐽 = (𝒖 𝑲(𝜌1𝑒 , 𝜌2𝑒 )𝒖 + 𝝀 (𝑲(𝜌1𝑒 , 𝜌2𝑒 )𝒖 − 𝒇)) (9)2𝜎𝑖𝑗,𝑗 = 0(1)

𝜎𝑖𝑗 = 𝐷𝑖𝑗𝑘𝑙 𝜖𝑘𝑙(2)1

𝜖𝑖𝑗 = (𝑢𝑖,𝑗 + 𝑢𝑗,𝑖 )2

Meanwhile, Eqs. (10) and (11) are obtained from thestationary condition of the Lagrange function. Eqs.
(12) and (13) can be obtained from Eqs. (10) and (11),respectively, and the corresponding equations for the
displacement and adjoint variable vectors (Eq. (14)are derived from Eq. (12) and (13), respectively. The
relationship derived is known as the self-adjoint
relationship.

(3)

Here, 𝜎𝑖𝑗 , 𝜖𝑖𝑗 , 𝐷𝑖𝑗𝑘𝑙 , and 𝑢𝑖 represent thecomponents of stress and strain tensors, the elasticity
coefficient tensor, and the displacement component,
respectively. Applying the finite element Galerkinprocedure to discretize Eqs. (1)–(3), the finite
element equation shown in Eq. (4) is obtained4). Here,𝒖𝑒 and 𝒇𝑒 are the displacement and external force
vectors, respectively; 𝑲𝑒 is the stiffness matrix, as
expressed in Eq. (5); 𝑩𝑒 and 𝑫𝑒 denote the B matrix
and elasticity coefficient matrix, respectively.
𝑲𝑒 𝒖𝑒 = 𝒇𝑒(4)

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

(5)

𝑒

𝜕𝐽∗
= 𝑲𝒖 − 𝒇 = 𝟎
𝜕𝝀

(10)

𝜕𝐽∗
= 𝑲𝒖 + 𝝀𝑇 𝑲
𝜕𝒖

= 𝒇 + 𝑲𝑇 𝝀 = 𝟎
𝑲𝒖 = 𝒇𝑲𝑇 𝝀 = −𝒇𝒖 = −𝝀

(11)
(12)
(13)
(14)

Here, 𝜌1𝑒 and 𝜌2𝑒 are introduced as the design
variables, and the Young’s modulus 𝐸 is represented
by Eq. (15), where 𝐸0 , 𝐸1 , 𝐸2 , and 𝑝 indicate thelower limit of the Young’s modulus, the Young’s
modulus for materials A and B, and the penalty
parameter, respectively.

Assembling 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), where Γ1 and Γ2 indicate the Dirichlet andNeumann boundaries, respectively, and the hat
symbol indicates the specified value.

𝐸(𝜌1 , 𝜌2 )𝑲𝒖 = 𝒇̂ on Γ1
𝒖=𝒖{
𝒇 = 𝒇̂ on Γ2

(6)

𝑝

𝑝

𝑝

𝑝

= ((1 − 𝜌2 )𝐸1 + 𝐸2 𝜌2 ) 𝜌1 + 𝐸0 (1 − 𝜌1 )
𝑝 𝑝

𝑝

(7)

= (𝐸2 − 𝐸1 )𝜌1 𝜌2 + (𝐸1 − 𝐸0 )𝜌1 + 𝐸0 (15)

To solve the maximally stiff structural problem, theperformance function is defined as shown in Eq. (8).

These design variables indicate the ratio of thematerial density to the standard density; 𝜌1𝑒 and 𝜌2𝑒𝜌11

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

are expressed as(8)𝜌𝜌1𝑒 = (𝜌A𝑒 ) and 𝜌2𝑒 = (𝜌Be ) ,A0

B0

where 𝜌Ae and 𝜌Be represent the material densities
for materials A and B of the two-phase material,
respectively; and 𝜌A0 and 𝜌B0 indicate the values ofthe standard density for materials A and B,
respectively. The gradient of the Lagrange function
with respect to 𝜌1𝑒 and 𝜌2𝑒 is expressed by Eqs. (16)and (17), respectively. These equations are used to
update the design variables.

The problem is to determine the optimal topology tominimize the performance function 𝐽 . Considering
the assembled finite element equation as theconstraint condition of the performance function, theperformance function extended by the adjoint
variable is obtained as shown in Eq. (9), which is2

Density-based topology optimization for maximally stiff structure problem using Two-phase materials

𝜕𝐽∗
𝜕𝑲
𝜕𝑲
= 𝝀𝑇𝑒 𝑒 𝒖𝑒 = −𝒖𝑇𝑒 𝑒 𝒖𝑒
𝜕𝜌1𝑒
𝜕𝜌1𝑒
𝜕𝜌1𝑒

(16)

𝜕𝐽∗
𝜕𝑲
𝜕𝑲
= 𝝀𝑇𝑒 𝑒 𝒖𝑒 = −𝒖𝑇𝑒 𝑒 𝒖𝑒
𝜕𝜌2𝑒
𝜕𝜌2𝑒
𝜕𝜌2𝑒

(17)

optimality criteria method (see Eqs. (24) and (25)),
where 𝑘 and 𝜂 indicate the number of iterations andthe weighting factor, respectively.

𝜌1𝑒 (𝑘+1) = 𝜌1𝑒 (𝑘) (If the gradients in Eqs. (16) and (17) are directly
employed, then a checkerboard-like structure is
obtained as the optimized result5). Therefore, thefiltering method, that is, the sensitivity filter, isgenerally applied to control the complexity of thestructure (Eqs. (18) and (19)). The smoothed gradient

𝜕𝐽∗

(18)

𝜕𝐽∗

2

𝑤(𝑥𝑗 , 𝑦𝑗 ) = 𝑅 − √(𝑥𝑗 − 𝑥𝑒 ) + (𝑦𝑗 − 𝑦𝑒 )2

𝑉𝑡𝑜𝑡𝑎𝑙

∗ (𝑘)
̅̅̅
𝜕𝐽
)
𝜕𝜌2𝑒

(

−𝛬2 (𝑘) (

(𝑘)
𝜕𝑉
)
𝜕𝜌2𝑒

(24)
𝜂

)

(25)

①Input ofcomputational conditions(20)②Finite element
analysis③Computation of
performance function④Judgement of
convergence

Yes⑦Processing
to separate
into individualmaterials

No⑤Computation of,and application of sensitivity filter

element, the total volume, and the average of design
variables 𝜌1𝑒 and 𝜌2𝑒 at the 0th iteration, respectively.
In two-dimensional problems, unit thickness is
assumed.

𝑣𝑒 𝜌2 𝑒̅̅̅̅̅(0)𝑉2 = ∑𝑚𝑥− 𝜌2𝑒 ≤ 0
𝑒=1

)

Start

̅̅̅̅̅
̅̅̅̅̅
(0)
(0)
𝑉𝑡𝑜𝑡𝑎𝑙 , 𝜌1𝑒 , and 𝜌2𝑒 indicate the volume of each

𝑡𝑜𝑡𝑎𝑙

(𝑘)

(19)

In addition, volume constraint conditions areincluded in the Lagrange function. The modified
Lagrange function is expressed by Eq. (21), where 𝛬1and 𝛬2 indicate the adjoint variables, and 𝑉1 and 𝑉2
are the volume constraint conditions. 𝑉1 and 𝑉2 areexpressed by Eqs. (22) and (23), respectively; 𝑣𝑒 ,𝐿 = 𝐽∗ + 𝛬1 𝑉1 + 𝛬2 𝑉2
𝑣𝑒 𝜌1 𝑒̅̅̅̅̅(0)𝑉1 = ∑𝑚𝑥− 𝜌1𝑒 ≤ 0
𝑒=1 𝑉

𝜕𝑉
)
𝜕𝜌1𝑒

The computational flow is illustrated in Fig. 1. In
Eqs. (24) and (25), the adjoint variables 𝛬1 and 𝛬2are calculated using the bisection method. However,at the final iteration, the computed adjoint variables
𝛬1 and 𝛬2 are directly used, and the intermediatematerial is obtained at the junction points of materials
A and B. Therefore, at the 7th step shown in Fig. 1,
the material separation process is computed.
Definitions and examples of the material separation
are shown in Fig. 2. The values 𝑙1 and 𝑙2 arecomputed using the bisection method as well.

weighting coefficient shown in Eq. (20). The
smoothing effect is adjusted to the filtering radius 𝑅,
and 𝑚𝑥 indicates the number of elements.

∑𝑚𝑥
𝑒=1 𝑤(𝑥𝑗 ,𝑦𝑗 )𝜕𝜌 𝜌𝑗
𝜕𝐽̅∗
𝑗
=
𝜕𝜌2𝑒
𝜌2𝑒 ∑𝑚𝑥
𝑤(𝑥
,𝑦
𝑗 𝑗)
𝑒=1

−𝛬1 (𝑘) (𝜌2𝑒 (𝑘+1) = 𝜌2𝑒 (𝑘) (

𝜕𝐽̅∗
𝜕𝐽̅∗
values are calculated from 𝜕𝜌1𝑒 and 𝜕𝜌2𝑒 using the∑𝑚𝑥
𝑒=1 𝑤(𝑥𝑗 ,𝑦𝑗 )𝜕𝜌 𝜌𝑗
𝜕𝐽̅∗
𝑗
=
𝜕𝜌1𝑒
𝜌1𝑒 ∑𝑚𝑥
𝑒=1 𝑤(𝑥𝑗 ,𝑦𝑗 )

𝜂

∗ (𝑘)
̅̅̅
𝜕𝐽
)
𝜕𝜌1𝑒

(

⑥Update of design variables,

Fig.1 Computational flow.

(21)
(22)
(23)

Design variables 𝜌1𝑒 and 𝜌2𝑒 are updated using the3EndTakahiko KURAHASHI, Sakura ICHIKAWA, Masayuki KISHIDA, Mizuki IKARASHI
and Tetsuro IYAMA
materials A and B, respectively.
Definition
The numerical results for each case are shown inMaterial A (Absence)Figs. 4 and 5. Fig. 4 shows the variation in the
Material A (Presence)
performance function for each case, and Fig. 5 shows
Material B (Absence)
a comparison of the optimized topologies with andMaterial B (Presence)
without complete separation of design variables 𝜌1and 𝜌2 . In Fig. 5, the value of strain energy for eachFind adjoint variables Λ1(=l1) and Λ2(= l2)case is shown in the lower row; the result shown on
so as to satisfy the volume constraint condition
the left is that without the application of theunder the above definition
separation process to materials A and B, and theAfter the above computation
result shown in the right is that with the application.
The results show that the value of strain energy in thecase with separation is lower than that without
separation.

2023/1/5

Fig.3 Computational model.
Table 1 Computational conditions.
Total number of elements

2400

Total number of nodes

2501

Mesh size [m]■ Material A): E=E1
■ Material B): E=E2
■ Void: E=E0

1

Material A : Young’s modulus () [Pa]

Material B : Young’s modulus () [Pa]

Poisson’s ration ν

Fig.2 Definition and example of material separation.

3. Numerical experiments
In this study, numerical experiments were
performed on the MBB model shown in Fig. 3. The
computational conditions are presented in Table 1.
The value of the strain energy was compared between
the results obtained with and without processing toclassify them into individual materials. This
investigation was performed by changing the volume
reduction rate of material B: Case 1, 60%; Case 2,
70%; Case 3, 80%; and Case 4, 90%. The volume
reduction rate of material A was fixed at 50%. Design
variables 𝜌1 and 𝜌2 are the design parameters for

1
2
0.3

Design domain [m2]40×60Convergence criterion ε

1.0×10-5

Filtering radius R [m]1.5

Penalty parameter p

3

Weighting factor η

0.75

Move limit of non-dimensional density

0.05

Fig.4 Variation in performance function.
(X-axis: Number of iterations; Y-axis: Strain energy)4Density-based topology optimization for maximally stiff structure problem using Two-phase materials
Case1 (Volume reduction rate of material B

Distribution of(%)

Distribution of(%))

Complete separation ofmaterials A and B

11.543[J]Distribution of(%)

10.994[J]Case3 (Volume reduction rate of material B

Distribution of(%)

Case2 (Volume reduction rate of material B

Distribution of(%)

Complete separation ofmaterials A and B

12.449[J])11.799[J]Case4 (Volume reduction rate of material B

Complete separation ofmaterials A and B

13.946[J]Distribution of(%))

Distribution of(%)

13.175[J]Distribution of(%))

Complete separation ofmaterials A and B.

15.819[J]14.884[J]Fig. 5 Comparison of optimized topologies with and without complete separation of design variables 𝜌1 and 𝜌2 .

4. Conclusions>, (accessed on 4 January, 2023). (in Japanese)2)

In this study, numerical experiments involving
maximally stiff structure problems using two-phasematerials were performed via density-based topology
optimization, and the strain energy was compared forcases with and without separation for two-phasematerials. The MBB model was employed in thecomputational model. Results showed that the strainenergy of the optimized result with separation for thetwo-phase material was lower than that without
separation.
Notably, we previously developed a new update
method for the design parameter6). In the future, we
intend to apply the method to the present study.

JFE

Steel

Corporationfrom<https://www.jfe-

steel.co.jp/release/2019/01/190109.html >, (accessed on 4
January, 2023). (in Japanese)
3)

Y. Wada, T. Yamamoto, K. Yuge, A. Nakamoto, I. Kizaki
and M. Ueno : Multiphase Topology Optimization under
Various Mix Conditions, Transactions of the Society ofAutomotive Engineers of Japan, Vol.49, No.2, pp.329-334,
2018. (In Japanese.)4)

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

5)

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

M. Kishida and T. Kurahashi : Proposal of a modified
optimality criteria method for topology optimization
analysis in 3-dimensional dynamic oscillation problems,

Acknowledgements: The computations were
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.

International

NikkeiBusiness

Publications,

Inc.

from

for

Numerical

Methods

inEngineering, Vol.123, Issue3, pp.866-896, 2022.
(Received July 19, 2023)

References
1)

Journal<https://xtech.nikkei.com/atcl/nxt/column/18/01516/00004/

5