## 2 Materials and Methods  
The BESO method employs discrete variables (0/1) where only a small subset of elements require updates during iterations, and material properties are represented using just two types of elements (solid/void), resulting in low memory consumption. In contrast, the SIMP method utilizes continuous density variables, necessitating thousands of material phases to represent graded densities, which significantly increases memory overhead. For 3D models, BESO avoids the storage and computation of full-density matrices inherent to SIMP by selectively retaining elements based on sensitivity thresholds. This approach mitigates the computational time escalation observed in SIMP for million-element models. Meanwhile, the MMC method relies on predefined parametric components, limiting its ability to freely generate complex internal cavities or curved boundaries. Currently, MMC remains relatively constrained in describing intricate topologies, requiring further methodological advancements. The BESO method employs finite element analysis (FEA) to define the relative density of elements as design variables. For discrete 0-1 optimization problems, a design variable value of 1 indicates a solid element, while a minimal value \(x_{\text{min}}\) (e.g., 0.001) represents a void element. This study adopts the soft-kill strategy [12] where material removal is achieved by transitioning an element's relative density from 1 to a near-zero value \(x_{\text{min}}\) (e.g., 0.001) rather than complete deletion. The mathematical formulation of BESO for compliance minimization under a prescribed volume fraction constraint is expressed as Formula (1):  
\[ \begin{split}& find\;\mathbf{x}=\left(x_{1},x_{2},\ldots x_{n_{t}} \right)\\ &\min\;\mathbf{C}=\tfrac{1}{2}\mathbf{U}^{T}\mathbf{K}\mathbf{U}=\tfrac{1}{2}x_{i} ^{p}\sum\limits_{i=1}^{n}u_{i}^{T}k_{0}u_{i}\\ & s.t.\;\mathbf{F}=\mathbf{K}\mathbf{U}\\ & V_{rep}-\sum\limits_{i=1}^{n}v_{i}x_{i}=0\\ & x_{i}=x_{\text{min}}\;or\;1\end{split}  \tag{1}\]  
where \(x_{i}\) denotes the element density variable, \(x_{i}=1\) represents solid elements and \(x_{i}=0.001\) denotes void elements, \(\mathbf{C}\) refers to the mean compliance, \(u_{i}\) is the element displacement vector, \(\mathbf{F}\) is the global force matrix, \(\mathbf{U}\) is the displacement matrix, \(\mathbf{K}\) is the global stiffness matrix, \(k_{0}\) is the element stiffness matrix, \(V_{rep}\) is the target volume fraction, and \(p\) is the penalization factor. The sensitivity analysis of the objective function is performed using the adjoint method [11], expressed as Formula (2):  
\[\alpha_{i}=-\frac{1}{p}\frac{\partial\mathbf{C}}{\partial x_{i}}= \begin{cases} \frac{1}{2}u_{i}^{T}K_{0}u_{i}&\text{when }x_{i}=1\\ \frac{p^{p-1}_{\text{min}}}{2}u_{i}^{T}K_{0}u_{i}&\text{when }x_{i}=x_{\text{min}} \end{cases}  \tag{2}\]  
where \(\alpha_{i}\) denotes the sensitivity value of the \(i^{th}\) element, quantifying its contribution to the objective function. The volume fraction at the current iteration step is computed via Formula (3) to establish the threshold for element addition/removal:  
\[V^{(l)}=\max\Bigl{\{}(1-c_{cr})V^{(l-1)},V_{rep}\Bigr{\}} \tag{3}\]where \(l\) is the number of iterations, \(c_{cr}\) is the evolutionary rate, and \(V^{(l)}\) is the Volume fraction to retain in each iteration. To mitigate checkerboard patterns and mesh dependency in topology optimization, a mesh-independent filtering scheme is implemented to regularize element sensitivity values, as expressed in Formula (4):  
\[\alpha_{e}=\frac{\sum\limits_{j=1}^{N}w_{ej}\alpha_{j}}{\sum\limits_{j=1}^{N}w_ {ej}} \tag{4}\]  
\[w_{ej}=\max\{0,r_{\min}-\Delta(e,j)\} \tag{5}\]  
where \(w_{ej}\) refers to the linear weighting factor, \(\Delta(e,j)\) is the distance between the centroids of element \(e\) and element \(j\), \(r_{\min}\) is the filter radius, and \(N\) is the total number of elements within the filter radius.