1. Introduction

1.1. Density-Based Topology Optimization

The density-based approach starts from a design domain Ω that is discretized using nodes and elements. This design domain is the maximum space (line, surface or volume) where the final design could exist, and a density value is assigned to each element of the discretized design domain, hence it is parametrized by this density. The density value is not binary, so it can vary from 0 (void) to 1 (solid), and this value is applied by using an interpolation function [8,9]. The use of only black and white elements easily leads to an ill-posed result [10]. Many solutions have been proposed [11,12,13], but none are totally satisfactory. An interesting proposal was to include not only binary values for the density of each element [9,14], those values being chosen automatically following a power rule (Zhou & Rozvany, 1991): the interpolation scheme, that will be explained later in this manuscript.

The used density-based approach is extensively described in several studies [6,16], and follows the simplified flowchart showed in Error! Reference source not found., where every step is numbered for future references.

Figure 1. Simplified flowchart of the density-based topology optimization process.

Figure 1. Simplified flowchart of the density-based topology optimization process.

In this study, special attention is given to levels 4 and 5 of the flowchart in Error! Reference source not found., as these processes involve the parameters under study. Specifically, for the interpolation in level 4, Solid Isotropic Material with Penalization method (SIMP) is utilized; additionally, a density filter is employed as the filtering technique in level 5.

1.2. Interpolation Scheme, the SIMP

The SIMP is based on the relation of the stress state of every element with its relative density. The formulation is included in Equation (1):

$${\mathrm{E}}_{\mathrm{i}}=\mathrm{E}\left({\mathrm{x}}_{\mathrm{i}}\right)={\mathrm{E}}_{\mathrm{m}\mathrm{i}\mathrm{n}}+{\mathrm{x}}_{\mathrm{i}}^{\mathrm{p}}({\mathrm{E}}_0-{\mathrm{E}}_{\mathrm{m}\mathrm{i}\mathrm{n}})$$

(1)

Being ${\mathrm{E}}_{\mathrm{i}}$ the relative Young Modulus of each element (calculated using FEM and representing the stress state of the element; note that it is also dependant on its degrees of freedom), ${\mathrm{E}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ as the almost zero elastic modulus of void material for avoiding mathematical incompatibility issues, ${\mathrm{E}}_0$ as the elastic modulus of the original material, ${\mathrm{x}}_{\mathrm{i}}$ as the relative density of the element $\mathrm{i}$, and $\mathrm{p}$ as the penalization factor. This penalization factor is one of the parameters to use in the sensitivity analysis of this manuscript. According to the Equation (1), it is evident that a higher $\mathrm{p}$ will result in a more binary outcome (Jiang, 2017), corresponding to the graph of a power rule function presented in Error! Reference source not found..

Figure 2. Graphical representation of the SIMP interpolation power rule for penalization values between 1 and 100.

Figure 2. Graphical representation of the SIMP interpolation power rule for penalization values between 1 and 100.

SIMP represents a robust, highly configurable, and computationally efficient scheme. It does not require homogenization and is adaptive to various design conditions and constraints. Due to these advantages, SIMP is chosen as the initial interpolation scheme for the research conducted in this work. However, using SIMP comes with challenges such as dealing with intermediate densities, mesh dependency, and potential non-convergence [17].

2. The Case Study

2.2. Case Study Description and Material

One aim of the study is to analyze parameters in a 3D structure, so the initial design space needs reliable dimensions in the x, y, and z directions. The typical load case involves minimizing compliance with a volume fraction constraint in a cantilever beam with a vertically applied distributed force. In this study, a point force (tip load) is used to emphasize the 3D nature of the case.

Considering that the element size is 1 in all directions, the initial design space consists of 80 elements in the x direction, 30 elements in the y direction, and 20 elements in the z direction. This choice helps generate fewer flat shapes in the optimization result. Error! Reference source not found. illustrates the case study with the applied load and the number of elements.

For this study, a volume fraction of 0.2 is applied, and a material with a Poisson Ratio (ν) of 0.3 is used. The value of ${\mathrm{E}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ in Equation 1 is set to 10⁻⁹ which is close enough to zero to avoid affecting calculations and prevent mathematical singularities; and ${\mathrm{E}}_0$ is set to 1, which is a typical used value in mentioned previous studies for researching and reducing calculation time. These properties are chosen for comparing the results to other studies [1,24] and having a direct application in the know and accesible MatLab^® codes [6,7].

Figure 3. Study Load Case and Design Domain.

Figure 3. Study Load Case and Design Domain.

The convergence is achieved when the maximum change in the density of any element between two consecutive iterations is less than 0.01, which is the recommended value based on the experience of previous studies [6].

2.3. Parameters for the Results Analysis

The sensitivity analysis aims to evaluate the machinability by exploring different filter radii and penalties to assist designers in their projects. The parameters extracted from the simulation results are intended to objectively demonstrate this characteristic. These parameters include the number of iterations, time per iteration, measure of non-discreteness, and machinability.

The number of iterations and time per iteration are directly obtained from MatLab^® and will be presented in the subsequent section. However, further explanations are needed for the measure of non-discreteness and machinability, which will be provided below.

2.3.1. Measure of Non-Discreteness

Discreteness refers to a global characteristic of the final design, indicating how binary the result is with distinct void (${\tilde{\mathrm{x}}}_{\mathrm{i}}$=0) and solid ${\tilde{\mathrm{x}}}_{\mathrm{i}}$=1) elements. This parameter is important because the representation of the result includes only elements with density ${\tilde{\mathrm{x}}}_{\mathrm{i}}$above a threshold (in this case, 0.5). Consequently, there are elements that do not appear in the final design but still contribute to compliance, and their significance increases as the discreteness of the result grows. Since the designer needs to interpret the shape displayed in the topology optimization (TO) result, and one objective of the study is to obtain easily interpretable components, it is essential to control and minimize the discreteness.

To analyze the characteristic of discreteness, Sigmund introduced a tool called the measure of non-discreteness (${\mathrm{M}}_{\mathrm{n}\mathrm{d}}$), which is widely utilized in the investigation of this matter [28]. The measure of non-discreteness is represented by Equation 6:

$${\mathrm{M}}_{\mathrm{n}\mathrm{d}}=\frac{{\sum }_{\mathrm{e}=1}^{\mathrm{N}}4{\overline{\mathrm{x}}}_{\mathrm{e}}(1-{\overline{\mathrm{x}}}_{\mathrm{e}})}{\mathrm{N}}\times 100\mathrm{\%}$$

(6)

In the provided context, $\mathrm{N}$ refers to the ID of each element in the design domain. The parameter ${\mathrm{M}}_{\mathrm{n}\mathrm{d}}$ represents the measure of non-discreteness. It should be noted that the lower the value of ${\mathrm{M}}_{\mathrm{n}\mathrm{d}}$, the more binary the result of the analysis is, indicating a higher level of discreteness.

2.3.2. Machinability

Here, the machinability of the optimized shape is analyzed by assessing the accessibility of a hypothetical tool to each frontier element. A frontier element refers to the last element with a density higher than the chosen threshold. This analysis involves an element-by-element evaluation.

To initiate the analysis, the set of elements exceeding the selected density threshold is divided into distinct groups:

• The external layer (EL): Comprising all the elements located on the boundaries of the design domain. These elements are always accessible because they do not have any solid interface with the exterior.
• The core layers (CL): Representing the internal solid region, consisting of all the elements forming the shape of the part excluding the frontier elements.
• The frontier layer (FL): Consisting of the elements above the density threshold and located between the solid and void phases.

The frontier layer (FL) can be determined by considering that it includes all the elements exceeding the threshold, but not belonging to the $EL$ or $CL$ groups, as stated in Equation (7):

$${\sum }_{\mathrm{n}=1}^{\mathrm{n}}\mathrm{N}={\sum }_{\mathrm{e}=1}^{\mathrm{e}}\mathrm{E}\mathrm{L}+{\sum }_{\mathrm{c}=1}^{\mathrm{c}}\mathrm{C}\mathrm{L}+{\sum }_{\mathrm{f}=1}^{\mathrm{f}}\mathrm{F}\mathrm{L}$$

(7)

The $\mathrm{F}\mathrm{L}$ is divided into three distinct groups in the subsequent step. These groups are based on the accessibility of the elements for a tool and their machinability. The groups are as follows:

• Elements directly accessible for a tool: These are the elements within the FL that can be readily machined using a tool.
• Elements not directly accessible but machinable by machining a neighboring element: These elements are not directly accessible to a tool but can still be manufactured through the machining of a neighboring element.
• Non-machinable elements: This group consists of elements within the FL that are not machinable.

This division allows for a more detailed understanding of the machinability characteristics of the elements within the FL, providing insights into their accessibility for machining operations. For the FL elements in groups 1 and 2, the process for identifying them is as follows:

For an element to be considered directly accessible for a tool, it must satisfy a specific condition: at its free faces, there must be a row of elements with densities below the selected threshold that extends to the boundaries of the design domain. To illustrate this process, let’s consider the simplified 2D example in Error! Reference source not found., where the identification process is depicted. The blue elements represent those with densities exceeding the threshold, while the green element e_(x,y) is the element currently being analyzed.

Figure 4. Element directly accessible for a tool detection process. In blue, solid phase; in green, the element under analysis in coordinates (x,y).

Figure 4. Element directly accessible for a tool detection process. In blue, solid phase; in green, the element under analysis in coordinates (x,y).

To determine if an element satisfies th of being directly accessible for a tool in the y axis, the row and column leading to that element are analysed from the outside towards the inside in both directions (shaded elements in the central draw of Error! Reference source not found.). If at any point during the analysis an element is encountered with a density exceeding the threshold, the condition isn’t fulfilled for that specific element.

To be considered machinable through its neighbor, an element must meet a specific condition: a free row must extend to its free face. To determine this, the neighboring rows are analyzed from the outside towards the inside until reaching the free face of the element. If, during this analysis, none of the studied lines contains an element with a density above the threshold, then the condition is fulfilled, indicating that the element can be machined through its neighboring elements. Error! Reference source not found. illustrates this process in a similar manner to the explained process in Error! Reference source not found. but in the x axis direction.

Figure 5. Element accessible for a tool through a neighbor detection process. In blue, solid phase; in green, the element under analysis in coordinates (x,y).

Figure 5. Element accessible for a tool through a neighbor detection process. In blue, solid phase; in green, the element under analysis in coordinates (x,y).

The output parameter is a percentage of the machinable elements with respect to all the frontier elements. The MATLAB ^® code of this process and the rest of the machining filter (that is explained in subsequent epigraphs) are shown in Appendix A.

1. 

Sensitivity Analysis

The results of the TOs are presented in Tables 1 to 4, which can be found in the Appendix B.1, Appendix B.2, Appendix B.3 and Appendix B.4. Note that the relative density ${\mathrm{x}}_{\mathrm{i}}$ of each element is represented with a greyscale factor in Tables 2 to 4; where a darker element represents a value closer to 1. Only elements with a relative density ${\mathrm{x}}_{\mathrm{i}}>0.5$ are shown. Finally, only full parts are considered, meaning that partial shapes with poor density concentration areas are discarded as satisfactory outcomes.

• Numerical results for the analysis are shown in Table 1, being:
• Sim. ID, the identifier of the specific TO case.
• Penalization, p; and Filter Radius R.
• Iterations, the number of iterations required to achieve the convergence condition. Note that a maximum of 5000 iterations was set, so if the TO reaches to this value, convergence is not reached.
• It. Time, the average time the TO took in every iteration.
• Objective, the compliance of the structure, calculated as in Equation (8):

$$\mathrm{c}\left(\mathrm{x}\right)={\mathrm{U}}^{\mathrm{T}}\mathrm{K}\mathrm{U}=\sum _{\mathrm{e}=1}^{\mathrm{N}}{\mathrm{E}}_{\mathrm{e}}\left({\mathrm{x}}_{\mathrm{e}}\right){\mathrm{u}}_{\mathrm{e}}^{\mathrm{T}}{\mathrm{k}}_0{\mathrm{u}}_{\mathrm{e}}$$

(8)

where $\mathrm{K}$, $\mathrm{U}$ and $\mathrm{F}$ are the global stiffness matrix and the global displacement and forces vectors, respectively; $\mathrm{N}$ is the total number of elements, ${\mathrm{u}}_{\mathrm{e}}$ is the e element displacement vector, ${\mathrm{k}}_0$ is the element stiffness matrix and ${\mathrm{x}}_{\mathrm{e}}$ is the e element density. It is calculated in the step 3 of Error! Reference source not found., the FEM analysis. It is used as a reference value for comparing different cases.

• Non-Discreteness, the described parameter in Section 2.3.1.
• Machinability, the described parameter in Section 2.3.2.

Next, in Error! Reference source not found. and Error! Reference source not found., some remarkable results are displayed along with a corrected graph of the discrete results data. The curve fitting strategy employed is a two-term exponential curve for Figure 6a, 6b and Figure 7a; and a three-term polynomial curve for Figure 7a.

As the filter radius increases, it is expected that higher filter radius values lead to simpler shapes, as the density filter in Equation (2) affects a greater number of elements. However, when the value of the radius becomes very large, the results start to exhibit discontinuities, and the computational cost increases.

Nevertheless, the corrected data depicted in Error! Reference source not found. reveals an unexpected trend. It shows that as the filter radius grows, the topology optimization (TO) process requires fewer iterations and less time per iteration for computation. This finding contradicts the initial expectation, as a larger filter radius entails analyzing more elements in each iteration.

Error! Reference source not found. provides an analysis of the manufacturability of the TOs. In Figure 7a, it is evident that higher penalization values lead to slightly lower M_nd. This is an expected outcome since higher penalization corresponds to a more aggressive interpolation of the SIMP function defined by Equation (1) and Error! Reference source not found..

Regarding non-discreteness, the results exhibit regular and predictable behavior for filter radii ranging from 4.25% to 6.5% of the EMD (Sim. ID 19-72), and the values increase with the corresponding radii. Similarly, machinability also increases with the filter radius since, as previously mentioned, a higher radius corresponds to simpler shapes that are more easily machinable.

In Figure 7, it is observed that results become more inconsistent and exhibit greater dispersion and poor predictability when large filter radii are used. Notably, in the default density-based TO approach, it is intriguing to note that achieving the best values for discreteness and machinability simultaneously is not possible, as depicted in Figure 8. The corrected data in Figure 8 is represented using a two-term exponential curve.

Machinability filter in combination with a Heaviside step filter [16,29,30,31] may result in an improvement of the ${\mathrm{M}}_{\mathrm{n}\mathrm{d}}$-Machinability combination, as will be exposed in the next epigraph.

Lower filter radius and lower penalization values lead to better compliance outcomes. Additionally, results show lower dispersion with the smallest radii. These findings are expected because more complex shapes result in higher moment of inertia, leading to increased stiffness. This can clearly be seen in Figure 9, where the corrected data is represented using a two-term exponential line.

2. 

Additions to Filtering

The filter is positioned inside the optimization algorithm, that in this case is the Optimality Criteria (OC) method described by Liu and Tovar in their work [6], using the parameters described in previous studies [32,33,34,35]. The flowchart depicted in Error! Reference source not found. illustrates the position of the filters (machining filter and Heaviside step filter) using the numbering notation from Error! Reference source not found..

Figure 10. Flowchart of the TO method with the machining filter highlighted in pink for clarity in later explanations.

Figure 10. Flowchart of the TO method with the machining filter highlighted in pink for clarity in later explanations.

Three additions are done to the code: the machining filter, the Heaviside step filter, and the ellipsoid-shaped density filter. These implementations are next explained and positioned inside the OC routine.

2.6. Ellipsoid-Shaped Density Filter

The main issue when extrapoling results from the sensitivity analysis to other design spaces is the mesh dependency of the density-based TO [17]. As the shape of the optimized part is directly dependant on the density filter described in Equation (2). This is caused by the relation of the filter size and the design domain size: using the same filter size with different design domains leads to different resultant shapes.

To adress this characteristic, the original sphere-shaped filter is taken as an ellipsoid with equal semi-axes, so they can be individually varied with the three cartesian dimmensions of the initial design domain. As the used weighting factor is conic weights [36], the radius of the ellipsoid at the same latitude and longitude than the analyzed element, ${\mathrm{R}}_{\mathrm{i}\mathrm{j}\mathrm{k}}^{\mathrm{f}\mathrm{i}\mathrm{l}}$, must be obtained for calculating the weight of the neighbor elements.

For obtaining this radius, six initial parameters are needed: $\mathrm{a},\mathrm{b}$ and $\mathrm{c}$, which are the semi-axes dimmensions given by the user as three filter radii ${\mathrm{R}}_{\mathrm{x}}^{\mathrm{m}\mathrm{i}\mathrm{n}},{\mathrm{R}}_{\mathrm{y}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ and ${\mathrm{R}}_{\mathrm{z}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ respectively; and the relative coordinates of the element being analysed in the Kernel filter with respect to the base element: ${\mathrm{R}}_{\mathrm{i}},{\mathrm{R}}_{\mathrm{j}}$ and ${\mathrm{R}}_{\mathrm{k}}$.

The angles $\mathsf{\alpha }$ and $\mathsf{\omega }$ in Error! Reference source not found. are the geocentric latitude and longitude respectively, they can be obtained with trigonometry, being:

$$\begin{array}{c}\mathsf{\alpha }={\tan }^{-1}\left({\mathrm{R}}_{\mathrm{y}}/{\mathrm{R}}_{\mathrm{x}}\right)\\ \mathsf{\omega }={\tan }^{-1}\left({\mathrm{R}}_{\mathrm{z}}/{\mathrm{R}}_{\mathrm{x}}\right)\end{array}$$

(11)

The coordinates of the ellipsoid surface at the same latitude and longitude than the analysed element, ${\mathrm{R}}_{\mathrm{x}},{\mathrm{R}}_{\mathrm{y}}$ and ${\mathrm{R}}_{\mathrm{z}}$, can be geometrically calculated. Combining these values with the general formula of an ellipsoid, Equation (12) is obtained for the calculation of the desired radius:

$${\mathrm{R}}_{\mathrm{i}\mathrm{j}\mathrm{k}}^{\mathrm{f}\mathrm{i}\mathrm{l}}=\frac{\mathrm{a}\mathrm{b}\mathrm{c}}{\sqrt{{\mathrm{c}}^2{\cos }^2\left(\mathsf{\omega }\right)\left({\mathrm{b}}^2{\cos }^2\left(\mathsf{\alpha }\right)+{\mathrm{a}}^2\mathrm{s}\mathrm{e}{\mathrm{n}}^2\left(\mathsf{\alpha }\right)\right)+{\mathrm{a}}^2{\mathrm{b}}^2{\cos }^2\left(\mathsf{\alpha }\right)\mathrm{s}\mathrm{e}{\mathrm{n}}^2\left(\mathsf{\omega }\right)}}$$

(12)

The three implementations are presented and highlighted in green in the flowchart presented in Error! Reference source not found.. The flowchart represents the level 6 of the Error! Reference source not found. and the blue block of the Error! Reference source not found.. The machining filter is divided into the three phases described in Appendix A.1, Appendix A.2 and Appendix A.3; and both density filters are computed with the ellipsoid-shaped neighborhood.

Figure 13. OC routine with the machining filter highlighted in pink and the colored in green.

Figure 13. OC routine with the machining filter highlighted in pink and the colored in green.

3. Results

References

1. Shiye, B.; Jiejiang, Z. Topology Optimization Design of 3D Continuum Structure with Reserved Hole Based on Variable Density Method. Journal of Engineering Science and Technology Review 2016, 9, 121–128. [CrossRef]
2. Michell, A.G.M. LVIII. The Limits of Economy of Material in Frame-Structures. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 1904, 8, 589–597. [CrossRef]
3. Dorn, W.S.; Gomory, R.E.; Greenberg, H.J. Automatic Design of Optimal Structures. Journal de Mécanique 1964, 3.
4. Bendsøe, M.P.; Kikuchi, N. Generating Optimal Topologies in Structural Design Using a Homogenization Method. Comput Methods Appl Mech Eng 1988, 71, 197–224. [CrossRef]
5. Sigmund, O.; Maute, K. Topology Optimization Approaches. Structural and Multidisciplinary Optimization 2013, 48, 1031–1055. [CrossRef]
6. Liu, K.; Tovar, A. An Efficient 3D Topology Optimization Code Written in Matlab. Structural and Multidisciplinary Optimization 2014, 50, 1175–1196. [CrossRef]
7. Ferrari, F.; Sigmund, O. A New Generation 99 Line Matlab Code for Compliance Topology Optimization and Its Extension to 3D. Structural and Multidisciplinary Optimization 2020, 62, 2211–2228. [CrossRef]
8. Bendsøe, M.P. Optimal Shape Design as a Material Distribution Problem. Structural Optimization 1989, 1, 193–202. [CrossRef]
9. Mlejnek, H.P. Some Aspects of the Genesis of Structures. Structural Optimization 1992, 5, 64–69. [CrossRef]
10. Kohn, R. V.; Strang, G. Optimal Design and Relaxation of Variational Problems, II. Commun Pure Appl Math 1986, 39, 139–182. [CrossRef]
11. Haber, R.B.; Jog, C.S.; Bendsøe, M.P. A New Approach to Variable-Topology Shape Design Using a Constraint on Perimeter. Structural Optimization 1996, 11, 1–12. [CrossRef]
12. Jog, C.S. A Robust Dual Algorithm for Topology Design of Structures in Discrete Variables. Int J Numer Methods Eng 2001, 50, 1607–1618. [CrossRef]
13. Duysinx, P. Layout Optimization: A Mathematical Programming Approach; Lyngby, 1997;
14. Bendsøe, M.P. Optimal Shape Design as a Material Distribution Problem. Structural Optimization 1989, 1, 193–202. [CrossRef]
15. Zhou, M.; Rozvany, G.I.N. The COC Algorithm, Part II: Topological, Geometrical and Generalized Shape Optimization. Comput Methods Appl Mech Eng 1991, 89, 309–336. [CrossRef]
16. Andreassen, E.; Clausen, A.; Schevenels, M.; Lazarov, B.S.; Sigmund, O. Efficient Topology Optimization in MATLAB Using 88 Lines of Code. Structural and Multidisciplinary Optimization 2011, 43, 1–16. [CrossRef]
17. Tyflopoulos, E.; Flem, D.T.; Steinert, M.; Olsen, A. State of the Art of Generative Design and Topology Optimization and Potential Research Needs.; Linköping, Sweden, August 14 2018.
18. Bendsøe, M.P.; Sigmund, O. Topology Optimization: Theory, Methods, and Applications.; 2nd ed.; Springer Science & Business Media, 2003;
19. Bruns, T.E.; Tortorelli, D.A. Topology Optimization of Non-Linear Elastic Structures and Compliant Mechanisms. Comput Methods Appl Mech Eng 2001, 190, 3443–3459. [CrossRef]
20. Altair Engineering Inc. Altair OptiStruct Help Guide.
21. Butze, M.; Sert, E.; Öchsner, A. Development of a Topology-optimized Indoor Crane Trolley for Additive Manufacturing. Materwiss Werksttech 2022, 53, 526–535. [CrossRef]
22. Dassault Systèmes SolidWorks Help Page.
23. PTC; Tonny Abbey PTC Creo Blogs.
24. Bendsøe, M.P.; Sigmund, O. Material Interpolation Schemes in Topology Optimization. Archive of Applied Mechanics (Ingenieur Archiv) 1999, 69, 635–654. [CrossRef]
25. Jiang, L.; Wu, C.W. Topology Optimization of Energy Storage Flywheel. Structural and Multidisciplinary Optimization 2017, 55, 1917–1925. [CrossRef]
26. Lee, H.Y.; Zhu, M.; Guest, J.K. Topology Optimization Considering Multi-Axis Machining Constraints Using Projection Methods. Comput Methods Appl Mech Eng 2022, 390, 114464. [CrossRef]
27. Langelaar, M. Topology Optimization for Multi-Axis Machining. Comput Methods Appl Mech Eng 2019, 351, 226–252. [CrossRef]
28. Sigmund, O. Morphology-Based Black and White Filters for Topology Optimization. Structural and Multidisciplinary Optimization 2007, 33, 401–424. [CrossRef]
29. Schevenels, M.; Sigmund, O. On the Implementation and Effectiveness of Morphological Close-Open and Open-Close Filters for Topology Optimization. Structural and Multidisciplinary Optimization 2016, 54, 15–21. [CrossRef]
30. Pellens, J.; Lombaert, G.; Lazarov, B.; Schevenels, M. Combined Length Scale and Overhang Angle Control in Minimum Compliance Topology Optimization for Additive Manufacturing. Structural and Multidisciplinary Optimization 2019, 59, 2005–2022. [CrossRef]
31. Wang, F.; Jensen, J.; Sigmund, O. Robust Topology Optimization of Photonic Crystal Waveguides with Tailored Dispersion Properties. JOSA B 2011, 28, 387–397. [CrossRef]
32. Prager, W. OPTIMALITY CRITERIA IN STRUCTURAL DESIGN. Proceedings of the National Academy of Sciences 1968, 61, 794–796. [CrossRef]
33. Karush, W. Minima of Functions of Several Variables with Inequalities as Side Conditions.; 2014.
34. Bendsøe, M.P. Optimization of Structural Topology, Shape, and Material; 1st ed.; Springer Berlin Heidelberg: Berlin, Heidelberg, 1995; ISBN 978-3-662-03117-9.
35. Sigmund, O. A 99 Line Topology Optimization Code Written in Matlab. Structural and Multidisciplinary Optimization 2001, 21, 120–127. [CrossRef]
36. Svanberg, K.; Svärd, H. Density Filters for Topology Optimization Based on the Pythagorean Means. Structural and Multidisciplinary Optimization 2013, 48, 859–875. [CrossRef]
37. Morillas, A.V.; Alonso, J.M.; Caballero, A.B.; Sisamón, C.C.; Ceruti, A. Adaptive Variable Design Algorithm for Improving Topology Optimization in Additive Manufacturing. Preprints (Basel) 2024. [CrossRef]