*Corr. Author’s Address: Wuhan University of Technology, School of Mechanical and Electronic Engineering, China, wuchaohua@whut.edu.cn 159
Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169 Received for review: 2023-07-25
© 2024 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2023-11-12
DOI:10.5545/sv-jme.2023.739 Original Scientific Paper Accepted for publication: 2024-01-25
An Eigenfrequency-Constrained Topology Optimization Method 
with Design Variable Reduction
Liu, W.C. – Wu, C.H. – Chen, X.G.
Wenchang Liu – Chaohua Wu* – Xingan Chen
Wuhan University of Technology, School of Mechanical and Electronic Engineering, China
The dynamic response of structures heavily relies on eigenfrequency, so the optimization of eigenfrequency is valuable in various working 
conditions. The bi-directional evolutionary structural optimization (BESO) method has been widely applied due to its ability to eliminate 
grayscale elements. Based upon BESO, this paper introduces a topology optimization method that incorporates eigenfrequency constraints 
and reduces the number of design variables. In this method, the optimization objective was to minimize compliance. The Lagrange multiplier 
was used to introduce eigenfrequency constraints, allowing for coordinated control of compliance and eigenfrequency. To prevent oscillation 
during the optimization process, the sensitivity was normalized. Additionally, to achieve faster convergence, the variables were reduced after 
meeting volume constraints. The numerical examples demonstrate the effectiveness of this method in increasing the eigenfrequency of the 
structure and avoiding resonance.
Keywords: Eigenfrequency constraint, topology optimization, bi-directional evolutionary structural optimization, design variable 
reduction; Lagrange multiplier method
Highlights
• 	 The eigenfrequency constraint was introduced through the Lagrange multiplier method.
• 	 To obtain faster convergence, the variables were reduced after meeting volume constraints.
• 	 The first-order natural frequency was increased by 42 % and 26.7 % in 2D numerical examples and 3D numerical examples 
respectively.
0  INTRODUCTION
Topology optimization is an optimization algorithm 
for material distribution. Compared with the 
traditional optimization methods, sufficient freedom 
is the biggest advantage of topology optimization, 
which provides a reliable and convenient solution for 
developing high-performance structures and obtaining 
the best material layouts [1] to [3].
After more than decades of research and 
development, the practicality of topology optimization 
has been fully proven. At present, continuum topology 
optimization methods mainly include Solid Isotropic 
Material with Penalization Method (SIMP), Level Eet 
Method (LSM), Evolutionary Structural Optimization 
(ESO), Bi-directional Evolutionary Structural 
Optimization (BESO), etc. [4] and [5].
The SIMP method was proposed by Bendsøe and 
Kikuchi [6], which based on the ideal of discretizing 
the design domain and relating the density and 
materials. This means that there is no material, or it 
is a solid material when the density value is 0 or 1. 
The variation of element density values between 0 
and 1 leads to the presenceof intermediate density 
elements, which is irrelevant in practical engineering. 
To address this issue, Rozvany et al. [7] proposed a 
density penalty scheme, which can update the element 
density towards 0 and 1 to obtain an approximate 0-1 
structure. Sigmund [8] proposed a sensitivity filtering 
method to eliminate the checkerboard pattern and 
mesh dependency, which makes SIMP more stable. 
Osher and Sethian [9] proposed the concept of a level 
set function, Sethian and Wiegmann [10] first to apply 
this method to topology optimization; it updates the 
structure through the continuous evolution of the level 
set function and obtains clear and smooth boundaries. 
The method has a slow convergence speed, and it is 
not easy to obtain hole structures. The ESO method 
was proposed by Xie and Steven [11], which based on 
the ideal of gradually removing less efficient materials 
until the material requirements are met; its update 
concept is simple and clear, completely different from 
traditional mathematical programming algorithms. 
However, due to the possibility of mistakenly 
deleting elements during the optimization process to 
obtain local optima, Huang and Xie [12] improved 
this method to BESO. In this method, the sensitivity 
information of each element needs to be calculated 
and then sorted. The threshold is determined based 
on the volume fraction of each step. The elements 
with sensitivity numbers greater than the threshold 
are retained as solid elements, while elements with 
sensitivity numbers less than the threshold are deleted. 
Even if an element becomes a void element by 
deletion, its sensitivity information is still preserved 
and can be reinstated as a solid element in subsequent 

Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169
160 Liu, W.C. – Wu, C.H. – Chen, X.G.
iterations, ensuring that the optimization result is the 
optimal solution.
Topology optimization has also been widely 
applied to avoid structural resonance. Munk et al. 
[13] studied the problem of frequency topology 
optimization under dynamic loads, and proposed 
a topology optimization method that can enhance 
the selected frequency and the reduce gap between 
frequencies. Li et al. [14] proposed a modified 
frequency band-constrained Heaviside function, 
which is continuously differentiable and beneficial 
for sensitivity analysis. The numerical examples 
demonstrate that this method can maximize the 
eigenfrequency of the structure. Du et al. [15] carried 
out research on fault safety topology optimization 
based on independent continuous mapping (ICM) and 
expanded it to the realm of frequency optimization. 
Kang, et al. [16] proposed a topology optimization 
method for large-scale frequency constraints. Li, 
et al. [17] proposed a topology optimization method 
for frequency optimization with periodic structures 
and reduced the amount of calculation in the process 
of frequency optimization by utilizing dimension 
reduction technology. Leader, et al. [18] considered 
both stress and frequency constraints and utilized 
the Jacobi-Davidson eigenvalue solving method to 
solve the natural frequency problem. Wang et al. [19] 
established a dynamic topology optimization model for 
long-span continuum, effectively improving the first-
order frequency. Xu et al. [20] proposed a frequency 
optimization problem with casting constraints, which 
can effectively obtain convergent solutions when 
dealing with frequency maximization problems. 
Su and Liu [21] studied the topology optimization 
of a coupled stress continuum to maximize the 
eigenfrequency, and also demonstrated the influence 
of eigen length on the results of eigenfrequency 
optimization. Ferrari et al. [22] proposed a frequency 
optimization method linked to a multi-mesh 
eigenvalue solver, greatly saving computational costs. 
Guan, et al. [23] proposed a multi-constraint topology 
optimization method with stress, displacement, 
and frequency constraints. Kim et al. [24] applied 
the topology optimization method to increase the 
frequency and reduce the noise of the steel wheel. 
Duan et al. [25] proposed a topology optimization 
method that can increase the frequency in a limited 
way while meeting the manufacturing constraints. Oh 
et al. [26] proposed a topology optimization method to 
maximize the operating frequency range of hyperbolic 
elastic meta-material and explained the mechanical 
knowledge of the model in detail. Vicente et al. [27] 
proposed a parallel topology optimization method for 
frequency optimization to find the optimal layout of 
materials from both macro- and micro-perspectives.
In addition, how to reduce the number of iteration 
steps and accelerate the convergence process are 
also a problem that need to be solved in topology 
optimization. Zheng et al. [28] introduced a freedom 
reduction mechanism in topology optimization, 
effectively accelerating the convergence process and 
saving the calculation cost. Jia et al. [29] combined 
ESO with LSM, which can reduce the number of 
iteration steps by automatically generating holes in the 
low-strain energy region near the node. Lian et al. [30] 
added a hierarchical mesh refinement algorithm into 
the moving morphable component (MMC) algorithm 
to improve convergence speed. Joo and Jang [31] 
proposed a deep neural network topology optimization 
algorithm, which can improve the convergence speed 
by obtaining the history of intermediate designs. Li 
and Zhang [32] used high noise and unbiased random 
gradients to update design variables and expedite the 
convergence process. Du et al. [33] shared a set of 
efficient topology optimization Matlab codes, which 
resulted in faster convergence speeds by removing 
the freedom not belonging to the transmission path in 
the finite element analysis. Yang et al. [34] proposed 
an adaptive step size strategy that multiplies the 
speeds of different nodes by different step sizes, 
which can accelerate convergence and also reduce 
mesh dependency. According to the characteristics 
of the BESO algorithm, Lin et al. [35] proposed a 
dynamic evolution strategy to accelerate convergence 
in topology optimization. Ren, et al. [36] used faster 
model reduction methods to enhance convergence 
speed.
In this paper, an eigenfrequency-constrained 
topology optimization method with design variable 
reduction is proposed, which can rapidly converge 
while increasing the eigenfrequency. Numerical 
examples demonstrate the effectiveness of this 
method.
1  TOPOLOGY OPTIMIZATION  
WITH EIGENFREQUENCY CONSTRAINTS
1.1  Problem Statement
When topology optimization is applied to structural 
design, volume is usually taken as the constraint, and 
the minimum compliance is taken as the optimization 
objective. This reflects the fact that stiffness is an 
extremely important objective in traditional structural 
design concepts. However, the topology optimization 
model will be multi-objective, multi-constraint and 

Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169
161 An Eigenfrequency-Constrained Topology Optimization Method with Design Variable Reduction
include other related conditions to meet engineering 
requirements and address complex working 
conditions.
In recent years, constraints other than volume 
introduced in topology optimization can be broadly 
divided into two categories. One type is related to 
manufacturing, including maximum and minimum 
size constraints, connectivity constraints, hole size 
constraints, hole number constraints, inclination 
angle constraints, and self-supporting constraints; 
The other type of constraint is functionality, such as 
displacement constraints, stress-strain constraints, 
fatigue constraints, and damage constraints.
As is well known, avoiding resonance is one of 
the important design objectives in structural design; 
it can be avoided by increasing the eigenfrequency of 
the structure. BESO is widely applied due to its simple 
concept and clear boundaries. The BESO topology 
optimization model with eigenfrequency constraints 
can be described mathematically as follows,
 
find X
X


 

xxxx
C
im jn
ij 11 12 13
12 12
1
2
,,,
,, ..., ;, ,..., ;
min

U UK U
KU F
KM u
T
nn
ij
n
Vq V
x
WY
x
s.t.
or





















2
0
0
1
min
 


, (1)
where the X is the design variable, x
ij
 is the ij
th
 
elemental density with a value of either 1 for solid or 
x
min
 (0.001 in this paper) for void, compliance C(X) 
is an objective function, K is the global stiffness 
matrix, U is the global displacement vector, F is the 
force vector, M is the global mass matrix, u
n
 is the 
eigenvector corresponding to ω
n
, V
0
 is the initial 
volume of structure, q is reserved volume ratio, V is 
the final structural volume, ω
n
 is n
th
 natural frequency,  
WY
 
is frequency constraint value.
The following Rayleigh quotient indicate the 
relationship between ω
n
 and u
n
, as follows (Eq. (2)),
 
n
n
T
n
n
T
n
2

uK u
uM u
. (2)
1.2  Material Interpolation Scheme
The material interpolation scheme applied in 
calculating compliance is expressed as follows, Eqs. 
(3) and (4):
 EE x
ij ij
penal
=
0
, (3)
 KK
ij ij
penal
x =
0
, (4)
where the E
ij
 is the ij
th
 elemental Young’s modulus, 
E
0
 is the Young’s modulus of the solid element, penal 
(penal = 3 in this paper) is a value used for the density 
penalty. When penal   ≤  2,  there  is  a  lar ge  amount  of 
porous material, and the optimized structure cannot 
be manufactured. When penal   ≥  3.5,  there  is  no 
significant change in the final topology result. When 
penal   ≥  4,  it  will  make  the  cal culation  very  slowly . 
Therefore, penal = 3 in this paper. K
ij
 is the stiffness 
matrix of the ij
th
 element, K
0
 is the stiffness matrix of 
the solid element.
To avoid local vibration modes during finite 
element analysis and the solution of frequencies, the 
material interpolation scheme is defined as follows, in 
Eqs. (5) and (6):
  xx
ij ij


0
,   (5)
 Ex
xx
x
xx
ij
penal
penal ij
penal
ij
penal












minm in
min
1
1
 
E
0
,   (6)
where ρ
0
 and ρ(x
ij
) denote respectively the material 
density of the solid and ij
th
.
1.3  Sensitivity Analysis
1.3.1  Lagrange Multiplier Method
In BESO topology optimization, constraints other 
than volume can be added with Lagrange multiplier 
method; the objective function is expressed as follows, 
Eq. (7):
 fC WY
n
      ,, 0 (7)
where λ is the Lagrangian multiplier.
In the BESO method, the Lagrange multiplier 
method has been widely applied to solve multi-
constraint problems. For instance, Huang and Xie 
[37] utilized this method to address displacement 
constraints, while Fan et al. [38] employed it to tackle 
stress constraints. The Lagrange multiplier method is 
suitable for obtaining optimal solutions under multiple 
constraints. It is easy to perform the sensitivity 
analysis by introducing a Lagrange multiplier to 
incorporate inequality constraints.
1.3.2  Sensitivity Number
In the BESO method, it is necessary to sort the 
sensitivity of each element and then update the 

Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169
162 Liu, W.C. – Wu, C.H. – Chen, X.G.
variables by determining a threshold based on volume 
constraints. According to Eq. (7), objective function 
sensitivity can be obtained as follows, Eq. (8):
 








f
x
C
xx
ee
n
e


. (8)
According to Eqs. (1), (3) and (4), compliance 
sensitivity can be obtained, as follows, Eq. (9):
 




C
x
penal x
e
e
penal
e
T
e
1
2
1
0
UKU , (9)
where U
e
 is the displacement vector of e
th
. According 
to  Eqs.  (5)  and  (6),  Eqs.  (10)  and  (1 1)  can  be  obtained 
by the derivation calculus, as follows:
 



M
M
x
ij
0
, (10)
 






K
K
x
x
x
penal x
ij
penal ij
penal
1
1
1
0
min
min
, (11)
where M
0
 is the mass matrix of the solid element. 
Frequency sensitivity can obtained based on Eqs. (10) 
and (11), as follows, Eq. (12):















n
en
n
T
penal e
penal
n
x
x
x
penal x
1
2
1
1
1
0
2
0
uK M
min
min
 
u
n
. (12)
According to Eqs. (8), (9) and (12), complete 
objective function sensitivity can obtained, as follows, 
Eq. (13):







f
x
penal x
x
x
pe
e
e
penal
e
T
e
n
n
T
penal
1
2
1
2
1
1
1
0
UKU
u 

min
min
n nalx
e
penal
nn








1
0
2
0
KM u  . (13)
Finally, the sensitivity number can be obtained 
based on the sensitivity analysis, as follows, Eq. (14):
 
e
e
penal
f
x



1
, (14)
where α
e
 is the sensitivity number of e
th
.
1.3.3  Variable Update Principle
The filtering scheme can be used to avoid 
checkerboard patterns and mesh-dependency, as 
follows, Eq. (15):
 


ij lk
el km in el k
e
t
m
l
n
k
li kj
maxr
,
,,
,
   





22
11
0 


e el kl k
m
l
n
k
el k
,
,
,

 













 11
  (15)
where  ∆
ij,lk
 is the distance between the centers of 
elements lk and ij; μ
e,lk
 is a weight factor, r
min
 is a filter 
radius, t is current iteration steps. It is effective to 
ensure a smoother optimization process and improve 
the stability of the optimization model by averaging 
three historical sensitivity number for averaging, as 
follows, Eq. (16):
 
e
t
e
t
e
t
e
t
 
 12
,   (16)
where λ is 0 when frequency constraints are met. The 
objective function aims to minimize compliances, 
which is equivalent to the original model. λ can be 
updated a value that satisfies the constraints can be 
updated until the constraint is satisfied.
The update method for the Lagrange multiplier is 
expressed as follows, Eqs. (17) and (18):
 ss s
tt 
 
1
05 .,
min
 (17)
 
t
t
t
s
s





1
1
1
1
, (18)
where s
t
 is a constant with a value range of s
min
 to 
1, s
min
 is a very small positive number; s
t+1
 and λ
t+1
 
can be updated to 1 and 0, respectively, when the 
constraint is met.
An appropriate Lagrange multiplier updating 
strategy is of great importance for achieving speed and 
accuracy. Lagrange multiplier updating strategies may 
vary in different constraint problems, which requires 
specific analysis according to the individual problems. 
It is particularly associated with the sensitivity to 
Lagrange multipliers and the nonlinearity of the 
optimization model. The oscillation is normal when 
using the Lagrange multiplier updating strategy 
optimization process. However, convergence becomes 
difficult when faced with numerous and large 
oscillations.
Therefore, the normalization strategy needs to be 
adopted to avoid oscillation, as follows, Eq. (19):
 


e
t e
tt
tt



min
maxm in
, (19)

Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169
163 An Eigenfrequency-Constrained Topology Optimization Method with Design Variable Reduction
where α
min
t
 and α
max
t
 is the minimum and maximum 
sensitivity value in the t
th
 iteration step, respectively.
During the optimization process, the variation 
pattern of volume is expressed, as follows, Eq. (20):
 VV VE R
tt
  

max, ,
1
1 (20)
where V
t
 is volume value in the t
th
 iteration step, ER is 
the volume evolution rate.
The convergence condition is the value of five 
relative changes in compliance less than 0.01, as 
follows, Eq. (21):
 
 






f
cc
c
t
t
t
t
t
t
t
t
t
4
5
9
4
00 1 ., (21)
where ∇f   is  the  value  of  5  relative  changes  in 
compliance, c
t
 is the compliance in the t
th
 iteration 
step.
2  2D NUMERICAL EXAMPLES
As shown in Fig. 1, the design domain is a 180:90 
rectangular region for a prescribed volume fraction 
of V   =  50  %.  The  beam  is  simply  supported  at  both 
ends and vertically loaded (P = 10 N) in the middle 
of its lower edge. The rectangular design domain is 
divided into 180×90 four-node plane stress elements. 
Young’s modulus E = 1 MPa, the volume evolution 
rate ER = 0.01, filter radius r
min
 is twice the length 
of the element side, Poisson’s ratio ν = 0.3 and mass 
density ρ = 0.001 kg/m³. The optimization objective 
is to minimize the compliance while satisfying the 
constraint on the first-order eigenfrequency. 
Fig. 1.  The design domain of 2D numerical example
The optimization results can be obtained by 
setting different optimization parameters, as shown in 
Fig. 2.
The optimal topology without any eigenfrequency 
constraint is shown in Fig. 2a for comparison. When  
ω
1
  is  constrained   to  be  150  rad/s,  155  rad/s,  170  rad/s, 
178 rad/s, 188 rad/s, the resulting topologies are 
shown  in  Fig.  2b  to  f.  Their  compliances  are  547.9578, 
681.4366,  812.0577,  710.2652,  1096.741 1.  Their  first-
order  eigenfrequency  are  150.8144  rad/s,  155.5822 
rad/s, 175.5030 rad/s, 184.1 129 rad/s, 188.7140 rad/s.
From Fig. 2b to f, it can be seen that the first-
order eigenfrequency increases gradually after the 
eigenfrequency constraint is introduced, which 
can meet the constraint conditions. When the first-
order eigenfrequency increases gradually, the 
compliance also increases gradually. It can be seen 
that the stiffness is sacrificed while satisfying the 
eigenfrequency constraint. 
Without introducing eigenfrequency constraints, 
the optimization process is shown in Fig. 3a, When 
the eigenfrequency constraint is WY = 178 rad/s, the 
optimization process is shown in Fig. 3b.
From Fig. 3a, it can be seen that the compliance  
C continuously increased without eigenfrequency 
constraints as the material is continuously removed. 
When the volume constraint is satisfied, the 
compliance C reaches its maximum value, and finally 
tend to be stable. The first-order eigenfrequency will 
initially increase and then decrease. During the entire 
optimization process, ω
1
 will fluctuate obviously, but 
it will eventually stabilize. C is stable growth, and 
there are no noticeable oscillations throughout the 
entire optimization process.
It can be seen from Fig. 3b that the topology 
optimization in the direction of satisfying the 
eigenfrequency constraint is carried out firstly when 
the eigenfrequency constraint is WY = 178 rad/s. 
Then the topology optimization is then performed 
in the direction of minimum compliance once the 
structure satisfies the eigenfrequency constraint. The 
entire optimization process involves a coordinated 
optimization of eigenfrequency constraints and 
minimum compliance. The Lagrange multiplier plays 
a coordinating role. In the continuous coordination, the 
local optimal solution satisfying the eigenfrequency 
constraint and the minimum compliance is finally 
obtained. During the mid-term stage of the 
optimization process, there will be a large oscillation. 
With the removal of the material,  fluctuates near 
WY = 178 rad/s and eventually satisfies the constraints. 
At this time, the compliance  is also converging.
In this example, ω
1
 = 132.1478 rad/s can be 
increased to ω
1
 = 188.7140 rad/s. The ω
1
 is increased 
by  42  %.  In  practical  enginee ring  applications,  the 
value  of 42 % eigenfrequency  increase  is undoubtedly 
huge, which can effectively avoid structural resonance. 
When the design domain, material parameters and 
constraint conditions vary, the effect of frequency 
enhancement will be significantly different. However, 

Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169
164 Liu, W.C. – Wu, C.H. – Chen, X.G.
this method is undoubtedly a solution and can provide 
a reference for related designs.
3  DESIGN VARIABLE REDUCTION
In topology optimization, the large number of variables 
necessitates extensive calculations and results in a 
high degree of freedom. In topology optimization, 
a) WY = 0 rad/s; ω
1
 =132.1478 rad/s; C = 521.2869 b) WY = 150 rad/s; ω
1
 =150.8144 rad/s; C = 547.9578
c) WY = 155 rad/s; ω
1
 =155.5822 rad/s; =681.4366 d) WY = 170 rad/s; ω
1
 =175.5030 rad/s; C = 812.0577
e) WY = 178 rad/s; ω
1
 =184.1129 rad/s; C = 710.2652 f) WY = 188 rad/s; ω
1
 =188.7140 rad/s; C = 1096.7411
Fig. 2.  The optimization results of 2D numerical examples
a)     b) 
Fig. 3.  The optimization process of 2D numerical examples

Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169
165 An Eigenfrequency-Constrained Topology Optimization Method with Design Variable Reduction
certain variables converge early and reach a stable 
state. Similarly, taking the simply supported beam in 
Section 2 as an example, when the eigenfrequency 
constraint is not introduced, the changes in some 
variables are shown in Fig. 4. Where the X(i,j) is the 
ij
th
 elemental density.
It can be seen from Fig. 4 that different variables 
converge at different times: some converge early, and 
some converge later. This paper defines that when 
the  value  of  the  continuous  5-step  variable  does 
not change (that is, when Eq. (22) is satisfied), the 
variable is a stable variable, and the variable that does 
not satisfy Eq. (22) is a free variable.
     xx xxxx
ij
t
ij
t
ij
t
ij
t
ij
t

 1234
55 or
min
. (22)
When the topology optimization satisfies the 
volume constraint, the number of design variables 
can be reduced. This means that the stable variables 
remain unchanged and no longer participate in the 
variable update. The algorithm flowchart is shown in 
Fig.  5.  For  comparison,  the  simply  supported  beam 
in Section 3 is also utilized as an example. When 
a)            b) 
c)            d) 
e)            f) 
Fig. 4.  The change of different variables

Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169
166 Liu, W.C. – Wu, C.H. – Chen, X.G.
for convergence in both methods is essentially the 
same. However, there is a significant reduction in 
the number of iterations after the design variable is 
decreased. It has been proven that reducing the design 
variables can effectively decrease the number of 
iterations and speed up convergence.
Table 1.  Iterations required for convergence
Eigenfrequen 
constraint [rad/s]
lterations without 
variable reduction
lterations with variable 
reduction
WY = 150 82 75
WY = 155 83 79
WY = 170 78 75
WY = 178 88 79
WY = 188 81 77
Fig. 7.  The comparison diagram of iterations
4  3D NUMERICAL EXAMPLES
In the 2D example, the effectiveness of the 
eigenfrequency constrained topology optimization 
method with design variable reduction has been 
fully proved. Next, a 3D example is used for simple 
verification.
As shown in Fig. 8, the design domain is a 
30 : 20 : 10 cube region in which the degree of freedom 
of the intermediate nodes on both sides is restricted, 
and a concentrated load of 1 N is applied to the 
midpoint of the bottom surface.
Fig. 8.  The design domain of 3D numerical examples
the eigenfrequency constraint is WY = 178 rad/s, the 
optimization  process  is  shown  in  Fig.  6  using  the 
algorithm depicted in Fig. 5.
Fig. 5.  The algorithm flow chart
Fig. 6.  The optimization process by using  
the algorithm shown in Fig. 5
It  can  be  seen  from  Fig.  6  that  when  the  design 
variable reduction mechanism is introduced, the entire 
optimization process can still obtain the results that 
meet the constraints, and the number of oscillations is 
significantly reduced.
This paper conducts a large number of numerical 
examples to compare the iterations required for 
convergence without the introduction of design 
variable reduction mechanism and with the 
introduction of design variable reduction mechanism. 
The comparison results are shown in Table 1 and Fig. 
7.
It can be observed from Table 1 and Fig. 7 that 
the change trend of the number of iterations required 

Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169
167 An Eigenfrequency-Constrained Topology Optimization Method with Design Variable Reduction
The prescribed volume fraction is V   =  15  %. 
The cube design domain is divided into 30×20×10 
four node plane stress elements. Young’s modulus  
E = 1 MPa, volume evolution rate ER = 0.01, filter 
radius r
min
  is  1.5  times  the  element  side  length, 
Poisson’s ratio ν  =  0.3 and mass density ρ = 1 kg/m
3
. 
The optimization objective is to minimize the 
compliance with the first-order eigenfrequency 
constrained. The number of iterations is abbreviated 
as it. The optimization results of the 3D numerical 
examples are shown in Fig. 9.
The optimal topology without any eigenfrequency 
constraint is shown in Fig. 9a for comparison. When 
ω
1
  is  constrained   to  be  285  rad/s,  290  rad/s,  295  rad/s, 
300 rad/s, 330 rad/s, the resulting topologies are 
shown in Figs. 9b to f. Their compliances are 4.0722, 
4.0208,  4.0060,  4.1454,  5.1 104.  Their  first-order 
eigenfrequency  are  299.461 1  rad/s,  303.2690  rad/s, 
306.1478  rad/s,  31 1.2603  rad/s,  339.2250  rad/s.  Their 
iterations are 132, 129, 129, 138, 139.
It can be seen from Fig. 9 that the first-order 
eigenfrequency can be improved by introducing 
eigenfrequency constraints. In this example, 
ω
1
  =  267.6767  rad/s  can  be  incre ased  to  ω
1
  =  339.2250 
rad/s. The ω
1
  is  increased  by  26.7  %.  The  number  of 
iterations is within 140 steps, and the convergence 
speed is faster.
5  CONCLUSIONS
In this paper, an eigenfrequency constrained topology 
optimization method with design variable reduction 
is proposed. Based on BESO, the eigenfrequency 
constraint is introduced using the Lagrange multiplier, 
and the topology optimization is performed with the 
objective of minimizing compliance. After satisfying 
the volume constraint, the design variable is reduced, 
which can significantly decrease the number of 
iterations and expedite convergence. The first-order 
eigenfrequency  can  be  increase d  by  42  %  and  26.7 
a) WY = 0 rad/s; ω
1
 = 267.6767 rad/s; C = 3.6387; it = 200 b) WY = 285 rad/s;  ω
1
 = 299.4611 rad/s;  C = 4.0722; it = 197
c) WY = 290 rad/s; ω
1
 = 303.2690 rad/s; C = 4.0208; it = 198 d) WY = 295 rad/s;  ω
1
 = 306.1478 rad/s; C = 4.0060; it = 200
e) WY = 300 rad/s; ω
1
 = 311.2603 rad/s; C = 4.1454; it = 201 f) WY = 330 rad/s; ω
1
 = 339.2250 rad/s; C = 5.1104; it = 203
Fig. 9.  The optimization results of 3D numerical examples

Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169
168 Liu, W.C. – Wu, C.H. – Chen, X.G.
%  in  2D  numerical  examples  and  in  3D  numerical 
examples,respectively.
6  ACKNOWLEDGMENTS
This work was supported by Fundamental Research 
Funds for the Central Universities (2019Ⅲ138cG).
7  REFERENCES
[1] Rosen, D.W. (2016). A review of synthesis methods for additive 
manufacturing. Virtual and Physical Prototyping, vol. 11, no. 4, 
p. 305-317, DOI:10.1080/17452759.2016.1240208.
[2] Liu, J., Gaynor, A.T., Chen, S., Kang, Z., Suresh, K., Takezawa, 
A., Li, L., Kato, J., Tang, J., Wang, C.C.L., Cheng, L., Liang, X., To, 
A.C. (2018). Current and future trends in topology optimization 
for additive manufacturing. Structural and Multidisciplinary 
Optimization, vol. 57, no. 6, p. 2457-2483, DOI:10.1007/
s00158-018-1994-3.
[3] Meng, L., Zhang, W., Quan, D., Shi, G., Tang, L., Hou, Y., Piotr, 
B., Zhu, J., Gao, T. (2019). From topology optimization design 
to additive manufacturing: Today’s success and tomorrow’s 
roadmap. Archives of ComputationalMethods in Engineering, 
vol. 27, p. 805-830, DOI:10.1007/s11831-019-09331-1.
[4] Dilgen, C.B., Dilgen, S.B., Aage, N., Jensen, J.S.. (2019). 
Topology optimization of acoustic mechanical interaction 
problems: a comparative review. Structural and 
Multidisciplinary Optimization, vol. 60, no. 2, p. 779-801, 
DOI:10.1007/s00158-019-02236-4.
[5] Deaton, J.D., Grandhi, R.V. (2014). A survey of structural and 
multidisciplinary continuum topology optimization: post 2000. 
Structural and Multidisciplinary Optimization, vol. 49, no. 1, p. 
1-38, DOI:10.1007/s00158-013-0956-z.
[6] Bendsøe, M.P., Kikuchi, N. (1988). Generating optimal 
topologies in structural design using a homogenization 
method. Computer Methods in Applied Mechanics and 
Engineering, vol. 71, no. 2, p. 197-224, DOI:10.1016/0045-
7825(88)90086-2.
[7] Rozvany, G.I., Zhou, M., Birker, T. (1992). Generalized shape 
optimization without homogenization. Structural Optimization, 
vol. 4, p. 250-252, DOI:10.1007/BF01742754.
[8] Sigmund, O. (2001). A 99 line topology optimization code 
written in Matlab. Structural and Multidisciplinary Optimization, 
vol. 21, p. 120-127, DOI:10.1007/s001580050176.
[9] Osher, S., Sethian, J.A. (1988). Fronts propagating with 
curvature-dependent speed: Algorithms based on Hamilton-
Jacobi formulations. Journal of Computational Physics, vol. 79, 
no. 1, p. 12-49, DOI:10.1016/0021-9991(88)90002-2.
[10] Sethian, J.A., Wiegmann, A. (2000). Structural boundary 
design via level set and immersed interface methods. Journal 
of Computational Physics, vol. 163, no. 2, p. 489-528, 
DOI:10.1006/jcph.2000.6581.
[11] Xie, Y.M., Steven, G.P. (1993). A simple evolutionary procedure 
for structural optimization. Computers & Structures, vol. 49, 
no. 5, p. 885-896, DOI:10.1016/0045-7949(93)90035-C.
[12] Huang, X., Xie, Y.M. (2007). Convergent and mesh-
independent solutions for the bi-directional evolutionary 
structural optimization method. Finite Elements in Analysis 
and Design, vol. 43, no. 14, p. 1039-1049, DOI:10.1016/j.
finel.2007.06.006.
[13] Munk, D.J., Vio, G.A., Steven, G.P. (2017). Novel moving 
isosurface threshold technique for optimization of structures 
under dynamic loading. AIAA Journal, vol. 55, no. 2, p. 638-
651, DOI:10.2514/1.J054692.
[14] Li, Q., Wu, Q., Liu, J., He, J., Liu, S. (2021), Topology 
optimization of vibrating structures with frequency band 
constraints. Structural and Multidisciplinary Optimization, vol. 
63, p. 1203-1218, DOI:10.1007/s00158-020-02753-7.
[15] Du, J.Z., Meng, F.W., Guo, Y.H., Sui, Y.K. (2020). Fail-
safe topology optimization of continuum structures with 
fundamental frequency constraints based on the ICM method, 
Acta Mechanica Sinica, vol. 36, p. 1065-1077, DOI:10.1007/
s10409-020-00988-7.
[16] Kang, Z., He, J., Shi, L., Miao, Z. (2020), A method using 
successive iteration of analysis and design for large-scale 
topology optimization considering eigenfrequencies. Computer 
Methods in Applied Mechanics and Engineering, no. 362, p. 
112847, DOI:10.1016/j.cma.2020.112847.
[17] Li, M., Cheng, Z., Jia, G., Shi, Z. (2019), Dimension reduction 
and surrogate based topology optimization of periodic 
structures. Composite Structures, no. 229, p. 111385, 
DOI:10.1016/j.compstruct.2019.111385.
[18] Leader, M.K., Chin, T.W., Kennedy, G.J. (2019). High-resolution 
topology optimization with stress and natural frequency 
constraints. AIAA Journal, vol. 57, no. 8, p. 3562-3578, 
DOI:10.2514/1.J057777.
[19] Wang, Y., Qin, D., Wang, R., Zhao, H. (2021). Dynamic topology 
optimization of long-span continuum structures. Shock and 
Vibration, no. 2021, p. 1-9, DOI:10.1155/2021/4421298.
[20] Xu, B., Han, Y.S., Zhao, L., Xie. Y.M. (2019). Topology 
optimization of continuum structures for natural frequencies 
considering casting constraints. Engineering Optimization, vol. 
51, no. 6, p. 941-960, DOI:10.1080/0305215X.2018.1506771.
[21] Su, W., Liu, S. (2016). Topology design for maximization of 
fundamental frequency of couple-stress continuum, Structural 
and Multidisciplinary Optimization, no. 53, p. 395-408, 
DOI:10.1007/s00158-015-1316-y.
[22] Ferrari, F., Lazarov, B.S., Sigmund. O. (2018). Eigenvalue 
topology optimization via efficient multilevel solution of the 
frequency response. International Journal for Numerical 
Methods in Engineering, vol, 115, no. 7, p. 872-892, 
DOI:10.1002/nme.5829.
[23] Guan, H., Chen, Y.J., Loo, Y.C., Xie, Y.M. (2003). Bridge 
topology optimisation with stress, displacement and frequency 
constraints. Computers & Structures, vol. 81, no. 3, p. 131-
145, DOI:10.1016/S0045-7949(02)00440-6.
[24] Kim, J., Kim, J.J., Jang, I.G. (2022). Integrated topology and 
shape optimization of the five-spoke steel wheel to improve 
the natural frequency. Structural and Multidisciplinary 
Optimization, vol. 65, no. 3, p. 78, DOI:10.1007/s00158-022-
03183-3.
[25] Duan, Z., Yan, J., Lee, I., Lund, E., Wang, G. (2019), Discrete 
material selection and structural topology optimization of 
composite frames for maximum fundamental frequency with 
manufacturing constraints. Structural and Multidisciplinary 

Strojniški vestnik - Journal of Mechanical Engineering 70(2024)3-4, 159-169
169 An Eigenfrequency-Constrained Topology Optimization Method with Design Variable Reduction
Optimization, no. 60, p. 1741-1758, DOI:10.1007/s00158-019-
02397-2.
[26] Oh, J.H., Ahn, Y.K., Kim, Y.Y. (2015). Maximization of operating 
frequency ranges of hyperbolic elastic metamaterials by 
topology optimization. Structural and Multidisciplinary 
Optimization, no. 52, p. 1023-1040, DOI:10.1007/s00158-
015-1288-y.
[27] Vicente, W.M., Zuo, Z.H., Pavanello. R., Calixto, T.K.L., Picelli., 
R., Xie, Y.M. (2016). Concurrent topology optimization for 
minimizing frequency responses of two-level hierarchical 
structures. Computer Methods in Applied Mechanics 
and Engineering, no. 301, p. 116-136, DOI:10.1016/j.
cma.2015.12.012.
[28] Zheng, W., Wang, Y., Zheng, Y., Da, D. (2020). Efficient topology 
optimization based on DOF reduction and convergence 
acceleration methods. Advances in Engineering Software, no. 
149, p. 102890, DOI:10.1016/j.advengsoft.2020.102890.
[29] Jia, H., Beom, H.G., Wang, Y., Lin, S., Liu, B. (2011). 
Evolutionary level set method for structural topology 
optimization. Computers & Structures, vol. 89, no. 5-6, p. 445-
454, DOI:10.1016/j.compstruc.2010.11.003.
[30] Lian, R., Jing, S., He, Z., Shi, Z., Song, G. (2020). An 
accelerating convergence rate method for moving morphable 
components. Mathematical Problems in Engineering, vol. 
2020, p. 1-15, DOI:10.1155/2020/2478292.
[31] Joo, Y., Yu, Y., Jang, I.G. (2021). Unit module-based 
convergence acceleration for topology optimization using the 
spatiotemporal deep neural network. IEEE Access, no. 9, p. 
149766-149779, DOI:10.1109/ACCESS.2021.3125014.
[32] Li, W., Zhang, X.S. (2021). Momentum-based accelerated 
mirror descent stochastic approximation for robust topology 
optimization under stochastic loads. International Journal for 
Numerical Methods in Engineering, vol. 122, no. 17, p. 4431-
4457, DOI:10.1002/nme.6672.
[33] Du, Z., Cui, T., Liu, C., Zhang, W., Guo, Y., Guo X. (2022). An 
efficient and easy-to-extend Matlab code of the Moving 
Morphable Component (MMC) method for three-dimensional 
topology optimization. Structural and Multidisciplinary 
Optimization, vol. 65, no. 5, p. 158, DOI:10.1007/s00158-022-
03239-4.
[34] Yang, C.D., Feng, J.H., Shen, Y.D. (2022). Step-size adaptive 
parametric level set method for structural topology 
optimization, Journal of Mechanical Science and Technology, 
no. 36, p. 5153-5164, DOI:10.1007/s12206-022-0928-6.
[35] Lin, H., Xu, A., Misra, A., Zhao, R. (2020). An ANSYS APDL code 
for topology optimization of structures with multi-constraints 
using the BESO method with dynamic evolution rate (DER-
BESO).  Structural and Multidisciplinary Optimization, no. 62, 
p. 2229-2254, DOI:10.1007/s00158-020-02588-2.
[36] Ren, C., Min, H., Ma, T., Wang F. (2020). Efficient structure 
crash topology optimization strategy using a model order 
reduction method combined with equivalent static loads, 
Proceedings of the Institution of Mechanical Engineers, Part 
D: Journal of Automobile Engineering, vol. 234, no. 7, p. 1897-
1911, DOI:10.1177/0954407019893841.
[37] Huang, X., Xie, Y.M. (2010). Evolutionary topology optimization 
of continuum structures with an additional displacement 
constraint. Structural and Multidisciplinary Optimization, vol. 
40, no. 1-6, p. 409-416, DOI:10.1007/s00158-009-0382-4.
[38] Fan, Z., Xia, L., Lai, W., Xia, Q., Shi, T. (2019). Evolutionary 
topology optimization of continuum structures with stress 
constraints. Structural and Multidisciplinary Optimization, no. 
59, p. 647-658, DOI:10.1007/s00158-018-2090-4.