*2.1. Topology Optimization Method*

The variable density topology optimization problem [32] can be expressed as follows:

$$\begin{aligned} \min\_{\rho\_{\varepsilon}} & \phi(\rho) = \mathcal{U}^{T} K \mathcal{U} = \sum\_{\varepsilon=1}^{N} \left(\rho\_{\varepsilon}\right)^{P} u\_{\varepsilon}^{T} k\_{0} u\_{\varepsilon} \\ \text{s.t.} & \begin{cases} K \mathcal{U} = F \\ \sum\limits\_{\varepsilon=1}^{N} \rho\_{\varepsilon} V\_{\varepsilon} \le V^{\*}, 0 \le \rho\_{\min} \le \rho\_{\varepsilon} \le 1 \end{cases} \end{aligned} \tag{1}$$

where *φ* means a sum of each element's compliance. *U*, *K* and *F* mean the global displacement, stiffness matrix and force vectors, respectively. *k*<sup>0</sup> and *ue* mean the element stiffness matrix and displacement vector. *ρ* means the design variable vector, viz element density vector. *N* means the element amount of design domain. *Ve* and *V*∗ mean the unit element volume and total volume for the design domain. To make sure of the numerically stable iteration, *ρ*min = 0.001 is chosen as the lower limitation of design variable. Additionally, a convolution-type filtering operation is used to filter the holes, viz.

$$\rho\_{\varepsilon} = \frac{\sum\_{i \in N\_{\varepsilon}} w(r\_i, r\_{\varepsilon}) v\_i x\_i}{\sum\_{i \in N\_{\varepsilon}} w(r\_i, r\_{\varepsilon}) v\_i} \tag{2}$$

where *Ne* = { *i*| *ri* − *re* ≤ *R*} means a neighborhood set within the radius, *R*. *ri* and *re* mean the spatial central coordinate of elements *i* and *e*, respectively. *w*(*ri*,*re*) = *R* − *ri* − *re* means the weighting function. *vi* is the volume of element *i*.
