On Non-Penalization SEMDOT Using Discrete Variable Sensitivities

This work proposes a non-penalization Smooth-Edged Material Distribution for Optimizing Topology (SEMDOT) algorithm, which is a typical elemental volume fraction-based topology optimization method, by adopting discrete variable sensitivities for solid, void, and assumed boundary elements instead of the continuous variable sensitivities used in the penalization one. In the proposed non-penalized SEMDOT algorithm, the material penalization scheme is eliminated. The efficiency, effectiveness, and general applicability of the proposed non-penalized algorithm are demonstrated in three case studies containing compliance minimization, compliant mechanism design, and heat conduction problems, as well as thorough comparisons with the penalized algorithm. In addition, the length scale control approach is used to solve the discontinuous boundary issue observed in thin and long structural features. The numerical results show that the convergency of the newly proposed non-penalization algorithm is stronger than the penalization algorithm, and improved results can be obtained by the non-penalized algorithm.


Introduction
Topology optimization is a light-weight design tool that can find the optimal material layout within a design domain subjected to given set of loads and boundary conditions [48]. One of the main issues associated with topology optimization is that the boundaries are not smooth because they are based on the edges of the finite element mesh. The challenge is to develop an efficient topology optimization method that can utilize a smooth boundary during the optimization.
Structural topology optimization can be traced back to the work by Cheng and Olhoff [8] and has been widely investigated since the pioneering paper regarding the homogenization method by Bendsøe and Kikuchi [5]. Some representative algorithms were developed, including the Solid Isotropic Material with Penalization (SIMP) [49], Evolutionary Structural Optimization (ESO) [60], Bi-directional Evolutionary Structural Optimization (BESO) [23], level-set method [57], Moving Morphable Component (MMC)-based method [20], and floating projection topology optimization algorithm [21,24]. These algorithms are widely used in solving different optimization problems, including large-scale problems [1,44], fluid and heat transfer problems [2], thermoelastic problems [19], multi-material problems [40], failure problems [39], and design for additive manufacturing [6,15,25]. In addition, articles regarding topology optimization algorithms with open access codes are comprehensively discussed and summarized by Wang et al. [55] for educational and academic purposes.
Despite their easy implementation and ability to remove materials freely across the design domain, element-based topology optimization algorithms will inevitably generate zigzag (discrete methods, such as BESO) or both zigzag and fuzzy boundaries (continuum methods, such as SIMP). Hence, post-processing smoothing methods are generally required to form smooth part surfaces after topology optimization to enable fabrication [29]. The post-processing methods will not only require extra effort, but it will increase the structural volume and mass when the maximum smooth part envelope is created from the topology optimization boundary [36,47]. Strategies that enable the formation of smooth topological boundaries in existing element-based algorithms have attracted a lot of research attention. Early works on smooth topological design were presented by Maute and Ramm [42], where a method that uses cubic or Bézier splines to generate smooth boundaries. Nguyen et al. [43] proposed a multi-resolution topology optimization algorithm that can form smooth boundaries via three distinct meshes (the displacement mesh, density mesh, and design variable mesh). Wang et al. [56] combined an adaptive mesh-adjustment method with SIMP to represent irregular boundaries using isoparametric elements. Using the BESO framework, Da et al. [9] developed an Evolutionary Topology Optimization (ETO) algorithm where smooth boundaries were generated based on the solid/void design of the grid points assigned to each element. This ETO method was extended by Liu et al. [35] to solve stressbased optimization problems. The effectiveness of the proposed algorithm was further verified by Chen et al. [7] and Li et al. [27] through implementing topology optimization of photonic crystals. Andreasen et al. [3] merged the level-set method into the density-based method where a cut element method is used to facilitate the generation of clear boundaries. Rodriguez et al. [46] incorporated Non-Uniform Rational Basis Spline (NURBS) hyper-surfaces into SIMP through which smooth topologies can be obtained. Very recently, Li et al. [28] proposed a boundary density evolution method without penalization where smooth boundaries are generated using the element refinement and level-set values based on the nodal strain energy. Qiu et al. [45] combined BESO and isogeometric analysis to obtain smooth topological designs by intersecting NURBS hyper-surfaces, describing the global sensitivity function, with suitable hyper-planes. Zhuang et al. [67] integrated a body-fitted triangular/tetrahedral mesh generation algorithm into BESO to match the dynamically updated boundaries exactly and form smooth topological boundaries.
Following the smooth-edged material distribution strategy presented in the ETO method and Huang's idea on smooth topological design [7,27], Fu et al. [11,12] proposed a new topology optimization tool named Smooth-Edged Material Distribution for Optimizing Topology (SEMDOT) based on the SIMP framework. SEMDOT performs at two levels: grid-points and elements. The grid-point values represent the solids, voids, and boundaries of the part(s) in the design envelope. A smooth boundary can be constructed from the grid-point data. There are a fixed number of grid-points that are associated with each element. The elements are used in the optimization analysis. It should be noted that SEMDOT is a relatively new topology optimization method in contrast to the pioneering algorithms of SIMP and BESO. The obvious advantages of SEMDOT are that it can be easily used with some existing methods developed based on SIMP and is able to form smooth boundaries in topology optimization without requiring any post-processing methods. The advantage of using smooth boundaries within the topology optimization method is that manufacturing and production constraints can be considered, implying that the SEMDOT method can be used for more than the initial design/prototyping concept stage. Fu et al. [10,13,14] integrated the additive manufacturing filter proposed by Langelaar [25] into SEMDOT to obtain 2D and 3D smooth self-supporting topologies. The algorithm mechanism of SEMDOT and its applications for additive manufacturing are systematically presented in the work by Fu [16].
In recent years, the SEMDOT method has broadened its range of applications to solve a wide range of optimization problems. Zhang et al. [65] adopted SEMDOT to solve natural convection heat transfer problems based on the reduced-order model. Yi et al. [63,64] used SEMDOT to generate smooth topologies for the research on improving energy performance in product designs for additive manufacturing. Gonçalves et al. [18] took advantage of SEMDOT to obtain smooth boundaries for the design of mechanical heterogeneous specimens by solving the compliant mechanism design problem.
The SEMDOT method is being continuously improved in the current version of SEMDOT [11] where the penalty coefficient of 1.5 is used instead of the "magic" coefficient of 3 used in SIMP. This is because the power-law model, a material penalization scheme, is applied to grid points in SEMDOT; however, SIMP uses the power-law model on elements to force the separation between solids and voids. SEMDOT pursues the solid/void separation at the grid point level based on both the power-law model and a Heaviside smooth function. Therefore, there is the potential to obtain solid/void grid points entirely based on using the Heaviside smooth function in SEMDOT. One large advantage of non-penalization topology optimization is that instead of using an artificial power-law relation, the physically meaningful relationship between the elemental value and its material properties can be implemented. Hence, the original optimization problem is not artificially changed to avoid imprecise solutions. The use of the material penalization scheme increases the nonlinearity and non-convexity of the objective functional or constraint function, which may cause the optimization algorithm to easily run into local optima [52]. The material penalization scheme may cause the overestimation of structural stiffness in density-based methods [59], and in terms of new optimization problems, the physical meaning of the material penalization model needs to be investigated and understood, because of the nonlinear scaling of the material and its effect on the multi-material structural design [21,26]. If the material penalization scheme can be removed from the algorithm, improved solutions can be expected, and multi-material optimization problems can be attempted. In addition, one parameter (the penalty coefficient) can be removed to reduce the algorithm complexity, and therefore, the discussion on the selection of a proper penalty parameter is not required in non-penalized algorithms. To date, no study has investigated the removal of the material penalization scheme from elemental volume fraction-based algorithms.
This motivates the presented work on developing the non-penalization SEMDOT method. The aim of this work is to further improve the capabilities of SEMDOT in structural performance and convergency for different topology optimization problems. This work will utilize discrete variable sensitivities for solid, void, and assumed boundary elements to replace the continuous variable sensitivities (based on the power-law model in SIMP) to realize the non-penalization strategy. Compared to some existing non-penalization strategies using specific filters [17,28], the non-penalization method here is proposed by slightly modifying the form of the sensitivity analysis, which enables SEMDOT to use more filters for various optimization purposes. As previous works regarding SEMDOT only gave the sensitivity analysis equations without detailed explanations [11,12], this work will present the details on the reason for using the assumed linear combination form for the sensitivity analysis. This work will take into account various representative topology optimization problems, including compliance minimization, compliant mechanism design, and heat conduction cases, to show the advantages of the non-penalized algorithm, and most of numerical cases presented here have not been discussed in previous works.

Algorithm Framework
The development of the SEMDOT algorithm takes full advantage of the benefits of the smooth representation in the ETO algorithm and element-based optimization in the SIMP algorithm. Despite being developed based on the ETO algorithm, SEMDOT is able to use both the optimality criteria method and Method of Moving Asymptote (MMA) as optimizers instead of the optimality criteria-based optimizer in the ETO, which enables SEMDOT to obtain the improved structural performance.
SEMDOT operates at two levels, grid-points and elements. The grid-point values represent the material distribution and the boundaries within the design volume. The elements are used to perform the optimization analysis, such as compliance minimization. SEMDOT generates a smooth topological boundary based on the solid/void design of the grid points that are independently assigned to each element [11], as illustrated in Fig. 1a. Unlike the SIMP algorithm that pursues the separation between solid and void elements, the basic idea of SEMDOT is to pursue the separation between solid and void at the grid-point level. Element properties are determined by whether they are closest to being solid, boundary, or void, as shown in Fig. 1b-d. This study focuses on the use of the structured finite element mesh to demonstrate the algorithmic mechanism of SEMDOT, and how the solid, void, and boundary elements are generated is based on the fixed-mesh finite element analysis.
In the penalized SEMDOT method [16], the power-law model in SIMP is applied to grid points, and hence, the Young's modulus of a grid point can be expressed by: where E e (ρ e,g ) is the function of the Young's modulus with respect to grid point densities, E 1 is the Young's modulus of the solid material, ρ e,g is the scaling value of the gth grid point assigned to the eth element, which is named the grid point density from here on (see Fig. 1) rather than the physical density, and ρ min is a small artificial value (e.g., 0.001). It is noted that ρ e,g is assigned 1 to represent a solid grid point or ρ min to represent a void grid point. In topology optimization, the exponent in the power-law function (i.e., p in Eq. (1)) is defined as the penalty coefficient.
In the new non-penalized SEMDOT method, the material penalization scheme is removed from the calculation of Young's moduli, and hence, Eq. (1) can be rewritten as: E e (ρ e,g ) = ρ e,g E 1 . (2) Because of using the SIMP framework at the element level, element-based variables are needed in SEMDOT to conduct finite element analysis. Elemental volume fractions that depend on grid point densities are defined as surrogate design variables such that where X e is the volume fraction of the eth element and N is the total number of grid points in each element. The function of Eq. (3) is to convert discrete variables (grid point densities ρ e,g ∈ {ρ min , 1}) to continuous variables (elemental volume fractions X e ∈ [ρ min , 1]). As shown in Fig. 1b, c, grid point densities are homogeneously distributed across solid (with a relative value of X e = 1) or void (with a relative value of X e = ρ min ) elements. Based on Eqs. (1) and (3), the Young's modulus of a solid or void element in penalized SEMDOT [16] is expressed by: Using the same form, the stiffness matrix of a solid or void element in penalized SEMDOT is: where K 1 e is the stiffness matrix of the solid element. Based on Eqs. (2) and (3), the Young's modulus of a solid or void element in non-penalized SEMDOT is The stiffness matrix of a solid or void element in non-penalized SEMDOT is therefore given by: Compared to the artificial relationship defined in Eq. (4), Eq. (5) establishes a physically meaningful relationship between the elemental volume fraction and its material property.
The difference between SIMP and SEMDOT is the intermediate elements. In SEMDOT, the boundary elements are artificially defined as the non-homogenized  Fig. 2b, the diagonal boundary across one element of X e = 0.5 is only for simplicity. During optimization, the arbitrary boundary associated with solid grid points can be formed across one element. Here, elemental stiffness matrices are estimated using a linear interpolation between the two phases of solid and void: where K 0 e is the stiffness matrix of the void element. The filter used for elemental volume fractions X e is whereX e is the filtered elemental volume fraction, N e is the neighborhood set of elements within the filter domain for the eth element that is a circle centered at the centroid of this element with a predefined filter radius r min , and ω el is a linear weight factor defined as where (e, l) is the center-to-center distance of the lth element within the filter domain to the eth element.
Nodal densities are obtained via a heuristic filter [23]: where ρ n is the density of the nth node, M is the total number of elements, and ω ne is the weight factor defined as where (n, e) is the distance between the nth node and the center of the eth element and ϒ min is the heuristic filter radius. It should be noted that the values of elemental volume fractions can be directly assigned to grid points without using the heuristic filter in Eq. (7). Given that the original aim of developing SEMDOT is to provide an easy-to-use and flexible design tool for the topology optimization community, it is recommended that the users keep the heuristic filter in Eq. (7) to yield a larger range of topological designs. Considering a four-node element, densities of grid points ρ(ζ, η) are calculated using the linear interpolation of nodal densities ρ n which has the form: where (ζ, η) is the local coordinate of the grid point, ρ γ n is the density for the γ th node of the element, and N γ (ζ, η) is an appropriate shape function which is expressed as [9]: where ζ γ and η γ are the local coordinates of the γ th node and ζ o and η o are the non-dimensional coordinates of the grid point o, as illustrated in Fig. 3.
The Heaviside smooth function is used to obtain solid/void grid points, which is where is a threshold value, (x, y) is the global coordinate of grid points, ρ(x, y) is the density of the grid point at (x, y), and β is a scaling parameter that controls the steepness and is updated by: where k is the current iteration number and is the evolution rate for β. In addition to the Heaviside smooth function in Eq. (8), other strategies that are capable of yielding solid/void grid points, for example, the Heaviside step function, can also be used. The thorough comparison between Heaviside smooth and step functions was presented in [11] where the advantages of using the Heaviside smooth function were comprehensively discussed. The smooth topological boundary is implicitly represented via a level-set function (x, y): where (x, y) is the level-set function for grid points. The threshold value in the level-set function (Eq. (9)) is determined in the Heaviside smooth function (Eq. (8)). It is noted that the value of is floating to meet the volume constraint during the optimization procedure.
For the next round of finite element analysis, filtered elemental volume fractions are updated by assembling grid points for their corresponding elements: where ρ new e,g is the updated density of the grid point obtained by the Heaviside smooth function.
In SEMDOT, there are two convergence criteria (i.e., the overall topological alteration and topological boundary error) used to terminate the optimization procedure. The overall topological alteration is defined as: where τ is the tolerance value for the overall topological alteration.
Solid element (14) Void element (4) Intermediate element on the boundary (5) Intermediate element not on the boundary (N v =2) Topological boundary

Fig. 4 Illustration of topological boundary error
The topological boundary error convergence criterion is defined as: where ε is the topological boundary error, N v is the number of intermediate elements that are not along boundaries, and is the tolerance value for the topological boundary error. The defined topological boundary error ε is schematically illustrated in Fig. 4 where a design domain with 25 elements is taken as an example. Equation (10) is an indicator used to assess the accuracy of the smooth topological boundary formed by the level-set function (Eq. (9)). Generally, increasing the number of elements within the design domain or grid points within each element is a simple way to improve the boundary accuracy, whereas the computational cost will be increased. The smooth-edged material distribution strategy proposed in SEMDOT could be wrongly applied to the final topology obtained by SIMP as the post-processing method. The formation of accurate smooth boundaries in SEMDOT requires a process of reducing the boundary deviation iteratively, which is assessed by the topological boundary error ε (Eq. (10)). A high boundary deviation related to Eq. (10) or discontinuous features could be caused if the smooth boundary strategy in SEMDOT is directly used as the post-processing method for SIMP. More details on the SEMDOT algorithm can be found in the works by Fu et al. [11,12].
with M design variables (elemental volume fractions X = (X 1 , . . . , X M )) subjected to a given iteration point X 0 is expressed by: where L e and U e are lower and upper moving asymptotes that are updated iteratively to mitigate oscillation or improve convergency, and the parameters r e and s e are defined as: Assuming compliance as the objective function which has the form of C(X e ) = f T u , for solid and void elements, the sensitivities with respect to X e in the penalization SEMDOT approach can be directly calculated based on Eq. (4): where C is the objective function and u e is the displacement vector of the eth element. Sensitivities of solid and void elements in SEMDOT (Eq. (11)) are calculated based on the idea of continuous variable sensitivities. An intermediate element (e.g., X e = 0.5) in SIMP is equivalent to an element with uniformly distributed intermediate grid points (i.e., ρ e,g = 0.5) in SEMDOT, as illustrated in Fig. 5 . Based on Eq. (4), the sensitivities of intermediate elements (ρ min < X e < 1) in SIMP are [4]: However, in SEMDOT, an intermediate element is artificially defined as a bimaterial boundary element by assembling solid (ρ e,g = 1) and void (ρ e,g = ρ min ) grid points (see Fig. 1d), and there are no intermediate grid points (0 < ρ e,g < 1) defined in boundary elements. That is, the sensitivity analysis in SIMP is based on the change of elements, whereas the sensitivity analysis in SEMDOT depends on the overall change of grid points (from 1 to 0 or from 0 to 1) across an element.
Based on the above discussion, Eq. (12) cannot be directly used to represent the sensitivities of boundary (intermediate) elements in SEMDOT. Here, physically meaningful sensitivities of boundary elements are calculated using a linear combination of Single elements of X e = 0.5 Intermediate grid point ρ e,g = 0.5 the sensitivities of void elements with a weighting factor of (1 − X e ) and sensitivities of solid elements with a weighting factor of X e , which has the form When the penalty coefficient p is set to 1 in Eqs. (11) and (13), we obtain ∂C(X e )/∂ X e = −u T e K 1 e u e for the sensitivities of all elements, which causes the difficulty in distinguishing solid, void, and boundary elements during the optimization process. Numerous intermediate elements will exist in topologies obtained by SEMDOT when using p = 1, and then, the optimization problem becomes the wellknown "variable-thickness-sheet" problem. The "variable-thickness-sheet" problem is a convex optimization problem, and therefore, the global optimum solution can be guaranteed [47]. However, SEMDOT does not establish the "variable-thickness-sheet" problem except for distinct structural designs.
To overcome the above issue caused when p = 1, the discrete variable sensitivity strategy presented in [30][31][32][33][34] is adopted through which sensitivities with respect to X e can be recalculated as: In Eq. (14), the effect of the penalty parameter p is completely removed, and hence, the non-penalization strategy is realized. The sensitivity analysis of void elements in Eq. (14) is a heuristic method. Discrete variable sensitivities, as shown in Eq. (14), are able to facilitate the distinction of solid and void elements based on their contribution to structural performance during optimization. The efficiency and effectiveness of the Fig. 6 Schematic illustration of sensitivities of boundary elements in non-penalized SEMDOT discrete variable sensitivity strategy are analytically and numerically validated in [34] by solving different topology optimization problems.
Unlike continuous variable sensitivities, discrete variable sensitivities are originally used in the sequential approximate integer programming algorithm, a discrete approach, and hence, there is no fixed form to calculate the sensitivities of intermediate elements. Based on discrete variable sensitivities in Eq. (14) and the linear combination in Eq. (13), the sensitivities of boundary elements in non-penalization SEMDOT can be written as: Equation (15) can be schematically illustrated by a plate with two different thicknesses (t 1 = 1 and t 2 = ρ min ) shown in Fig. 6, which is inspired by Sun et al. [53]. In Fig. 6, the whole plate ( ) is the combination of the plate with the thickness of t 1 = 1 ( 1 ) and the plate with the thickness of t 2 = ρ min ( 2 ).
Based on Fig. 6, sensitivities of boundary elements in non-penalized SEMDOT can be estimated by: where V 1 and V 2 are the volumes of domains 1 and 2 , respectively; V t is the volume of the plate with only one thickness of t 1 = 1; and A 1 , A 2 , and A are the areas of domains 1 , 2 , and , respectively. In Eq. (16), ∂C(X e )/∂ X e | X e =1 represents the sensitivity of the plate with only one thickness of t 1 = 1. Figure 6 can also be used to demonstrate the Young's moduli and stiffness matrices of boundary elements in non-penalized SEMDOT. By contrast, Eq. (13) in penalized SEMDOT cannot be schematically demonstrated by the plate shown in Fig. 6. Therefore, the newly developed non-penalized SEMDOT has more physical meaning than penalized SEMDOT.
If the sensitivities of boundary elements (ρ min < X e < 1) are calculated directly based on the assumed stiffness matrices of boundary elements related to Eq. (6) instead of using the linear combination, we obtain that As ρ min is a small value of 0.001, Eq. (17) is equivalent to ∂C(X e )/∂ X e = −u T e K 1 e u e . As previously discussed, Eq. (17) will cause the difficulty in forming proper topological layouts, and a large volume of intermediate elements will exist within the design domain. Therefore, it is not recommended to use Eq. (17) to calculate the sensitivities of boundary elements in SEMDOT.
Despite Eqs. (13) and (15) having the similar form, they are formulated based on two different ideas: continuous and discrete variable sensitivities. Equation (13) is calculated based on sensitivities of void and solid elements related to the power-law model in SIMP (Eq. (4)). By contrast, Eq. (15) is calculated based on discrete variable sensitivities of void and solid elements. It should be noted that material properties and sensitivity analysis of the boundary elements are not limited to the linear combination. Some nonlinear combinations or other strategies (for example, high-order sensitivity analysis) can also be considered to pursue more accurate sensitivities. Here, we only present the simplest form to reduce the algorithm complexity and to make it easy for potential users to understand the algorithm mechanism of SEMDOT.
It is admitted that the same form, = (1− X e ) Void + X e Solid , for the assumption of material properties and sensitivities of boundary elements causes a less rigorous mathematical expression. Mathematically, sensitivities are calculated based on the derivative of the function of the stiffness matrix with respect to the design variable. Compared to traditional element-based algorithms, the SEMDOT approach provides a new way to explore different topological designs. The effectiveness of the assumptions of sensitivity analysis in Eqs. (13) and (15) will be validated in the next section through assessing obtained topologies and objective function values.

Case Studies
This section will demonstrate the effectiveness and general applicability of the SEM-DOT method using discrete variable sensitivities in three case studies containing compliance minimization, compliant mechanism design, and heat conduction cases. Comparisons of SEMDOT with existing element-based topology optimization algorithms were comprehensively conducted by Fu et al. [11,12], and the effects of the parameters in SEMDOT on the results of optimization were thoroughly discussed by Fu et al. [11,12]. This study focuses on the comparison between non-penalization and penalization algorithms subjected to the same set of parameters. Unless otherwise stated, the parameters of SEMDOT in [11,12] and default MMA parameters are used. It is noted that the default penalty coefficient p is 1.5 in the penalized algorithm.

Compliance Minimization Case
The deep cantilever beam case depicted in Fig. 7 is applied to the compliance minimization or stiffness maximization optimization problem, which is stated as where f and u are global force and displacement vectors, respectively; K(X e ) is the global stiffness matrix; V e is the volume of the eth element; and V * is the target volume. In Fig. 7, L represents the length of the design domain, and F represents a unit vertical load. It is noted that L = 200, F = −1, r min = 3, and V * = 0.5 are used in this optimization case.
As shown in Fig. 8a, the optimization process of non-penalized SEMDOT related to Eq. (15) converges after 140 iterations, which is more efficient than that of penalized SEMDOT related to Eq. (13) (154 iterations). Although the structural compliance of penalized SEMDOT (60.8786) is slightly lower than that of non-penalized SEMDOT (60.9246), it cannot be concluded that non-penalized SEMDOT sacrifices the performance to improve the efficiency. For the fairness of comparison, the penalty parameter p should be removed from the calculation of the objective function value in penalized SEMDOT. The corrected objective function value of penalized SEMDOT C * (X e ) is 60.9651, which is slightly higher than that of non-penalized SEMDOT (60.9246). Two different smooth topological designs obtained by non-penalized and penalized SEMDOT are shown in Fig. 8b, c. The numerical investigation on penalization SEMDOT with p = 1 is conducted to explain the reason of proposing a separate non-penalization version of SEMDOT using the optimization case presented in Fig. 7. Figure 9a shows that the optimization process of penalized SEMDOT with p = 1 cannot converge after 300 iterations, with a structural compliance value of 55.9821 and a high topological boundary error value of 0.6836. In addition, the improper smooth topological design in Fig. 9b is caused by numerous intermediate elements shown in Fig. 9c. This demonstrates the high value of topological boundary error at the 300th iteration. Although the structural compliance of penalized SEMDOT with p = 1 is the lowest because of being the "variablethickness-sheet" problem, the smooth topological design is not obtained (see Fig. 9b). Therefore, the penalty parameter p should be greater than 1 ( p > 1) in penalized SEMDOT for convergency and reasonable topologies. The porous structure design problem is a suitable test case to evaluate the capability of SEMDOT in the aspect of forming smooth boundaries as numerous holes will be formed in this case [58]. Here, the porous structure constraint proposed by Wu et al. [58] is formulated by where is the targeted local volume fraction, q is the parameter that controls the accuracy of per-element constraints, andX e is defined as  where N e is the set of all surrounding elements with a centroid that is closer than a predefined radius R to the centroid of the eth element (X e ), which has the form: A finite element mesh of 400 × 200 (L = 400), q = 16, r min = 4, R = 28, V * = 0.5, andα = 0.65 are used in the porous structure design problem. Figure 10 a shows that there are obvious fluctuations at the first 150 iterations in both nonpenalized and penalized SEMDOT, and then, both the optimization processes reach steady state. The optimization process of non-penalized SEMDOT converges at the number of iterations of 337, which is 34.69% faster than that of penalized SEMDOT (516 iterations). The structural compliance of penalized SEMDOT is 65.4996, and its corrected value is 65.6269, which is 0.49% lower than that of non-penalized SEMDOT (65.9492). The difference of the objective function values between non-penalized and penalized SEMDOT can be ignored. The intact porous structure obtained by nonpenalized SEMDOT is shown in Fig. 10b. By contrast, penalized SEMDOT obtains a porous structure with broken thin features, as shown in Fig. 10c. It is noted that the broken thin features shown in Fig. 10c have slight contribution to the stiffness of the structure but will cause difficulties in manufacturing. For a more sophisticated optimization problem, non-penalized SEMDOT can lead to a more reasonable solution.

Compliant Mechanism Design Case
The second optimization case is the compliant mechanism design shown in Fig. 11 where a single-piece flexible structure transfers an input force to another point through the elastic deformation [22,50]. In Fig. 11, f in represents the input force; k in and k out represent input and output spring stiffnesses, respectively; and u out represents the output displacement. A mesh of 60 × 30 (L=60), f in = 1, k in = 0.1, k out = 0.1, r min = 1.5, and V * = 0.3 are applied to this case. In addition, the move limit in MMA is set to 0.5 to stabilize convergence. The objective of this optimization case is to achieve the maximum output displacement u out . The corresponding optimization problem is stated as: where L is a unit length vector with zeros at all degrees of freedom except at the output point where it is one,ũ is the dummy load displacement vector calculated by solving the adjoint problem Kũ = −L, and f in is the input force vector.
Sensitivities of the solid and void elements for compliant mechanism design in penalized SEMDOT are [50]: whereũ e is the dummy load displacement vector of the eth element. Sensitivities of the boundary elements in penalized SEMDOT are therefore estimated as: Sensitivities of the solid and void elements for compliant mechanism design in non-penalized SEMDOT are [32] ∂C(X e ) ∂ X e =  Sensitivities of the boundary elements in non-penalized SEMDOT are therefore approximated as: Figure 12a shows that the optimization process of non-penalized SEMDOT related to Eq. (19) converges at the output displacement of 0.9988 after 107 iterations, and the optimization process of penalized SEMDOT related to Eq. (18) converges at the output displacement of 1.0502 after 93 iterations. The corrected objective function value (u * out ) of penalized SEMDOT is 0.9944, which is lower than that of non-penalized SEMDOT (0.9988). In contrast to the performance, the convergency of penalized SEMDOT is slightly stronger than that of non-penalized SEMDOT in terms of this case. Smooth topologies obtained by non-penalized and penalized SEMDOT can be found in Fig. 12b, c where the main difference of these two structures is the size of the enclosed void area (hole).  In non-penalized SEMDOT, the Heaviside smooth function (Eq. (8)) is the only scheme used to facilitate the generation of solid/void elements, and therefore, the evolution rate in the Heaviside smooth function has effect on convergency, performance, and topological layout. Performance and convergency of non-penalized SEMDOT under six different evolution rates (i.e., = 0.125, 0.25, 0.75, 1, 1.5, and 2) are summarized in Table 1 . It is noted that the default evolution rate is 0.5 in SEMDOT. Although it is hard to identify the relation between the evolution rate and output displacement, a relatively small value of the evolution rate can lead to a better solution, for example, the output displacement of 1.0031 at = 0.125 in Table 1. As outlined in Table 1, increasing the evolution rate is able to accelerate the convergence process. The resulting smooth topologies under different evolution rates are shown in Fig. 13where no sharp or slender hinges are observed.
Generally, reducing the output spring stiffness k out will increase the difficulty in solving the force inverter design problem [41]. A more challenging case with k in = 1 and k out = 0.001 is involved in the comparison between non-penalized and penalized SEMDOT. In addition, a mesh of 100×50 (L=100), f in = 1, r min = 2.5, and V * = 0.3 are adopted. Figure 14a shows that the optimization process of non-penalized SEM-DOT converges at the output displacement of 1.4569 after 222 iterations, and penalized SEMDOT converges at the output displacement of 1.0148 after 583 iterations. In this challenging case, non-penalized SEMDOT performs much better than penalized SEMDOT in both convergence and performance. More importantly, non-penalized  SEMDOT yields a correct smooth topology (see Fig. 14b), whereas penalized SEM-DOT results in an improper topology (see Fig. 14c). Compared to penalized SEMDOT, non-penalized SEMDOT can deal with more difficult optimization problems.
To further challenge non-penalized SEMDOT, a case with k in = 1 and k out = 1 × 10 −9 is taken into account. When using k out = 1 × 10 −9 , a fine mesh that is capable of increasing the accuracy of the sensitivity analysis related to Eq. (19) is recommended in non-penalized SEMDOT. In this case, a mesh of 400 × 200 (L=400), f in = 1, r min = 5, and V * = 0.3 are used. Figure 15 a shows that the optimization process terminates at the number of iterations of 658 with the output displacement of 2.2113, and the topological boundary error reaches 0 after 157 iterations. The successfully optimized topology is shown in Fig. 15b.

Heat Conduction Case
The last case is the heat conduction problem, which is a challenge to the smooth boundary formation in SEMDOT due to the tree-like material distribution with thin and long features [38,61]. The typical heat conduction optimization problem is depicted in Fig. 16 where the square design domain is homogeneously heated, and a unique  heat sink with the temperature of zero (T = 0) is placed in the middle of the left edge. Two material phases with isotropic conductivities of 1 and 0.001 are taken into account here. An outer move limit on the maximum design variable alternation at each iteration is set to 0.01 to stabilize the optimization procedure of this 2D heat conduction case. In this case, a mesh of 120 × 120 (L=120), r min = 2, and V * = 0.5 are adopted. The goal of the considered case is to minimize the thermal compliance, which is equivalent to minimizing the average temperature. The optimization problem is mathematically described as: where ϕ is the thermal compliance, U is the nodal temperature vector, and P is the global thermal load vector which has the form of P =KU whereK is the thermal conductivity matrix.
Sensitivities of the solid and void elements for the heat conduction problem in penalized SEMDOT are [37]: where U e is the temperature of the eth element andK 1 e is the thermal conductivity matrix of the solid element.
Sensitivities of the boundary elements in penalized SEMDOT are expressed by: Sensitivities of the solid and void elements for the heat conduction problem in non-penalized SEMDOT are [62]: Sensitivities of the boundary elements in non-penalized SEMDOT are calculated as: Figure 17a shows that the optimization process of non-penalized SEMDOT related to Eq. (21) converges at the thermal compliance of 386.4578 after 375 iterations, and the optimized smooth topology is shown in Fig. 17b. The obvious oscillation is found in the optimization procedure of penalized SEMDOT related to Eq. (20), which causes the difficulty in convergence. Therefore, the optimization process of penalized SEMDOT is artificially terminated at the thermal compliance of 381.3003 after 1000 iterations. Due to the numerical instability, an improper topological design is obtained by penalized SEMDOT, as shown in Fig. 17c. The penalty coefficient p is increased from 1.5 to 3 to mitigate the aforementioned issue in penalized SEMDOT. Figure 18a shows that the optimization process of penalized SEMDOT with p = 3 converges at the thermal compliance of 392.6831 after 451 iterations, and the topological boundary error reaches a value close to zero. The corrected objective functional value of penalized SEMDOT with p = 3 is 411.8227, which is 6.56% higher than that of non-penalized SEMDOT (386.4578). Compared to the topology in Fig. 17c, an improved topology shown in Fig. 18b is obtained.
To demonstrate the generality of non-penalization SEMDOT in solving heat conduction problems, six different target volumes (i.e., V * =0.1, 0.2, 0.3, 0.4, 0.6, and 0.7) are considered, and corresponding thermal compliance and numbers of iterations   are listed in Table 2. As outlined in Table 2, increasing the target volume is a straightforward way to reduce the thermal compliance; however, more materials are required. It is hard to identify the relation between the total number of iterations and target volume. Optimized smooth topologies under distinct target volume fractions are shown in Fig. 19. The common issue of optimized topologies in Figs. 17b, c, 18b, and 19 is discontinuous boundaries at the tips of the obtained tree-like structures. The element-based topological design subjected to V * = 0.3 shown in Fig. 20 is taken as an exam- The discontinuous boundary issue in the heat conduction problem cannot be simply solved by increasing the number of elements or grid points. Although the two filters in SEMDOT can implicitly control the feature size, the discontinuous boundary issue still cannot be solved by increasing the two filter radii (i.e., r min and ϒ min ). Therefore, an effective approach is required to explicitly control the minimum feature size with the aim of overcoming the discontinuous boundary issue. The length-scale control approach proposed by Zhang et al. [66], which is originally developed for SIMP, is adopted here. In this method, the minimum length scale can be accurately measured based on the structural skeleton, which is identified by means of the image processing technique. The minimum length-scale constraint incorporated into SEMDOT is defined as: where X j is the elemental volume fraction of the jth element, X is the vector of X j , is a closed domain, j is the domain occupied by the jth element, SS( ) is the structural skeleton, d is the lower bound of the length scale, and B(X, d) is a closed ball centered at X with a diameter of d. Here, d is set to 4 element widths.
Again, the heat conduction problem subjected to V * = 0.3 is taken as an example to demonstrate the effectiveness of the selected length-scale control approach. Figure 21a shows that thermal compliance keeps increasing after 70 iterations, and hence, there is no steady state during optimization. In this case, the topological boundary error in Eq. (10) can be used to determine a proper topological design considering the minimum length-scale constraint. Figure 21b, c shows the structural skeleton and topology at the 100th iteration when ε is 0, and the resemblance between the topologies in Figs. 21c and 19c can be found. Converged skeleton and topology are shown in Fig. 21d, e, respectively. By contrast, the topology in Fig. 21e, which is obtained after 293 iterations, is far from that in Fig. 19c. Due to the use of the length scale constraint, 12.77% performance needs to be sacrificed to yield the topology in Fig. 21c, and 32.57% performance needs to sacrificed to obtain the topology in Fig. 21e. Solving the discontinuous boundary issue for SEMDOT is not the focus of this work, and hence, only one example (Fig. 21) is presented here.
Following the 2D heat conduction case, a 3D case sketched in Fig. 22 is taken into account to further show the advantage of non-penalized SEMDOT in solving 3D optimization problems. Here, a mesh of 60×60×60, r min = 4, V * = 0.1, and the move limit of 0.5 in MMA are utilized.
As discussed in the 2D case, increasing the value of the penalty parameter p can mitigate the oscillation of the optimization process in penalized SEMDOT. Therefore, the penalized SEMDOT method with p = 3 is used for comparison. Figure 23a shows that the optimization process of non-penalized SEMDOT converges at the thermal compliance of 824,899.956 after 800 iterations. By contrast, the optimization process of penalized SEMDOT with p = 3 is artificially terminated after 1500 iterations at the thermal compliance of 968,386.807. The corrected thermal compliance of the penalized SEMDOT with p = 3 is 1,382,263.89, which is 67.57% higher than that of non-penalized SEMDOT (824,899.956). The resulting smooth 3D topologies obtained by non-penalized and penalized SEMDOT are shown in Fig. 23b, c, respectively. Even though the length-scale control approach is not applied to the 3D heat conduction case, no discontinuous boundary issue is observed in Fig. 23b. This is because cutting an element in a 3D space can alleviate the feature breakage. However, discontinuous boundaries are observed in Fig. 23c. It is noted that the 3D smooth design in Fig. 23c is not the converged result. Based on the above discussion, it is concluded that, compared to penalized SEM-DOT, non-penalized SEMDOT can converge faster and yield better solutions for both 2D and 3D heat conduction problems.

Conclusions
The presented study developed a non-penalization SEMDOT algorithm without significant modifications based on discrete variable sensitivities. The efficiency, effectiveness, and general applicability of non-penalized SEMDOT are demonstrated via three representative case studies: compliance minimization, compliant mechanism design, and heat conduction problems. In addition, comparisons between penalized and non-penalized SEMDOT were thoroughly discussed. For the fairness of comparison, the objective function in penalized SEMDOT was corrected by removing the penalty parameter.
The proposed non-penalized SEMDOT method enables designers to directly obtain the smooth topologies without needing post-processing methods for engineering applications, and discussions on effect of the penalty coefficient p on final results are not required. Analytical and numerical discussions demonstrate that the assumed linear combination form for the sensitivity analysis in SEMDOT has physical meaning and it is effective. Thanks to discrete variable sensitivities, non-penalized SEMDOT can lead to improved solutions, stronger convergency, more reasonable topologies compared to penalized SEMDOT at least for the topology optimization problems presented in this work. Porous structure design and length-scale control cases prove that the existing methods developed based on SIMP can be easily incorporated into SEMDOT, which facilitates the use of SEMDOT for different optimization problems.
Future work will use the newly proposed non-penalized SEMDOT method to solve dynamic and wave propagation problems.

Declarations
Conflict of interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.