Level set band method: A combination of density-based and level set methods for the topology optimization of continuums

Abstract

The level set method (LSM), which is transplanted from the computer graphics field, has been successfully introduced into the structural topology optimization field for about two decades, but it still has not been widely applied to practical engineering problems as density-based methods do. One of the reasons is that it acts as a boundary evolution algorithm, which is not as flexible as density-based methods at controlling topology changes. In this study, a level set band method is proposed to overcome this drawback in handling topology changes in the level set framework. This scheme is proposed to improve the continuity of objective and constraint functions by incorporating one parameter, namely, level set band, to seamlessly combine LSM and density-based method to utilize their advantages. The proposed method demonstrates a flexible topology change by applying a certain size of the level set band and can converge to a clear boundary representation methodology. The method is easy to implement for improving existing LSMs and does not require the introduction of penalization or filtering factors that are prone to numerical issues. Several 2D and 3D numerical examples of compliance minimization problems are studied to illustrate the effects of the proposed method.

Key words： level set method; topology optimization; density-based method; level set band

Introduction

The successful topology optimization of continuums inspired from the study of the optimal thickness distribution of elastic plates [¹] has led to the successive development of several topology optimization models, such as the homogenization method [²], solid isotropic material/microstructure with penalization (SIMP) method [³–⁵], rational approximation of material property method [⁶], evolutionary structural optimization (ESO)/bidirectional ESO (BESO) [⁷–⁹], level set method (LSM) [¹⁰–¹⁴], independent-continuous mapping method [¹⁵], moving isosurface threshold method [¹⁶], stiffness spreading method [¹⁷,¹⁸], and moving morphable component/void method [¹⁹,²⁰]. The most mature method is SIMP, which has been successfully implemented in commercial software systems, such as OptiStruct, Tosca, and ANSYS. SIMP has been widely used worldwide, especially after Sigmund published the 99-line MATLAB code [²¹]. It provides an efficient way for new researchers to accept and start their work in topology optimization for continuums. ESO, which was first proposed by Xie and Steven [⁷], uses a similar idea to that of SIMP. The idea is to find an appropriate way of material distribution in the design domain. Later, Refs. [⁸,⁹,²²] proposed BESO, which solves the drawback of ESO that a material cannot be added back into a structure after it is deleted. ESO/BESO has a clear concept of updating a structure to achieve optimum design, and it is easy to understand and implement. It has also been widely studied and implemented in the topology optimization tool Ameba.

LSM [¹⁰–¹⁴,²³–³⁰], which is borrowed from computer graphics and image-processing fields, was first implemented in topology optimization around 2000 [¹¹,¹²]. It uses a different idea based on the method of front tracking to drive the boundary of a structure iteratively to obtain the optimum design. It has received wide attention since it was first introduced into the topology optimization field and became popular in a short time due to its advantages [¹³,¹⁴]. For example, it always provides clear structural boundaries or material interfaces, making it suitable for optimization problems related to geometric control, and it does not require the introduction of a penalty factor in contrast to SIMP; therefore, it is more stable when solving dynamic optimization problems [³¹]. However, LSM is developed from a method based on the boundary evolution concept. The iterative updating of a structure during optimization is always generated from boundaries rather than from the entire domain as in density-based methods. Consequently, it lacks nucleation capacity and involves strong initial design-dependent problems.

To overcome these drawbacks, we analyze the fundamental scheme of LSM to evaluate its difficulties in handling topology changes and compare the differences between the conventional LSM (CLSM) and the variational LSM, which is called zero LSM (ZLSM) in this paper. The comparison demonstrates that the latter is more recommended for topology optimization in practical applications because it has the ability of nucleation and less stringent requirements on meshes compared with the former. A level set band method is proposed to improve the capability of LSMs in handling topology changes. The method is easy to implement by involving only one parameter, the level set band, and does not require the introduction of penalization or filtering factors that are prone to numerical issues. This proposal may pave the way for the wide acceptance of level set-based topology optimization method in practical engineering applications.

The remainder of this paper is organized as follows. The basic concept of LSM is introduced and a comparison between CLSM and ZLSM is presented in Section 2. The level set band method, which is a new method with a variable level set band, is proposed in Section 3. The optimization model is introduced in Section 4 and the method is studied and evaluated with several numerical examples in Section 5. Finally, conclusions are provided in Section 6.

LSM for topology optimization

LSM was first proposed by Osher and Sethian [¹⁰] for interface tracking. It has been successfully implemented in computer graphics and image segmentation fields before introducing into the topology optimization field [²³,²⁴]. The concept of a 2D design, which is represented with a level set function, is illustrated in Fig. 1.

Fig.1 Basic concept representing a boundary \partial \Omega with a zero level set.

Full size|PPT slide

The basic idea of LSM is to define a one dimension higher surface, i.e., the level set function \mathrm{\Phi }, and use the zero level set of the function to represent the boundary of an object and control the evolution of the boundary by updating the level set surface, as illustrated in Fig. 1. The definitions of the boundary and the different parts of the domain are given as follows:

(1)

where \mathit{x}\in R^{\mathrm{2}} or R^{\mathrm{3}} denotes a point in the design domain D\subset R^{\mathrm{2}} or R^{\mathrm{3}}, \Omega and \partial \Omega are the solid part and boundary of the object, respectively, and t is the pseudo time to represent the updating iteration steps during the evolution of the level set function.

Conventional level set method

In CLSM, i.e., the earliest proposed LSM [²³,²⁴], the name “LSM” consists of a systematic methodology to trace the front of an interface in an implicit way in a fixed Euler grid mesh. This methodology includes a discrete way to solve the Hamilton–Jacobi (H–J) equation by implementing upwind scheme, ENO, or WENO [³²,³³], the re-initialization process to recreate a signed distance level set function [²⁴,³⁴], and the velocity extension approach [³⁵] to extend the boundary velocity to the entire evolution domain or a narrow band area around the boundary to alleviate the computation cost. This methodology can be considered a conventional way to implement LSM in related problems. The approaches have been summarized in two classical books [²³,²⁴]; LSM has been well developed in computer graphics, as demonstrated by the brilliant animations and its successful implementation in the film industry. However, the application of LSM in the field of topology optimization has remained at the academic level for a long time, and its potential has not been fully exploited due to several limitations. In the following part, the basic concept of CLSM is schematically discussed, and its inherent limitations are analyzed.

In CLSM, a level set function is updated by solving H–J equation given in Eq. (2), which is a partial differential equation with spatial derivative \nabla \mathrm{\Phi } and temporal derivative \partial \mathrm{\Phi }/\partial t:

(2)\dot{\mathrm{\Phi }}\equiv {\displaystyle \frac{\partial \mathrm{\Phi }}{\partial t}}-V_{\mathrm{n}}\left|\nabla \mathrm{\Phi }\right|=\mathrm{0},

where V_{\mathrm{n}} in topology optimization is usually obtained using the sensitivity analysis based on the theory of shape derivative by applying the finite element method (FEM) to solve the state and adjoint equations [¹⁴]. The “-” in Eq. (2) depends on the outside positive normal direction of V_{\mathrm{n}}, which is determined by the level set representation model given in Eq. (1). If \mathrm{\Phi }(\mathit{x},\mathrm{\hspace{0.17em}}t)\hspace{0.17em}>\hspace{0.17em}\mathrm{0} is defined as the outside part, then “-” should be changed to “\mathrm{+}”.

Although CLSM has been introduced into the topology optimization field for a long time, it is still not as well understood by most people as density-based methods. The updating scheme of CLSM should be demonstrated with a diagram, as shown in Fig. 2, where the level set function is updated with one small time step \mathrm{\Delta }t. \dot{\mathrm{\Phi }} in Eq. (2) consists of two parts, which are \left(\partial \mathrm{\Phi }/\partial t\right)\mathrm{\Delta }t and V_{\mathrm{n}}|\nabla \mathrm{\Phi }|\mathrm{\Delta }t. The two parts are illustrated in Fig. 2(a). The equation \dot{\mathrm{\Phi }}=\mathrm{0} is always satisfied in each step. Therefore, the two parts are cancelled each other out, and a point on the level set surface \mathrm{\Phi } of the nth time step pⁿ moves horizontally to point pⁿ⁺¹ in the next time step. Here, V_{\mathrm{n}} denotes the normal velocity on point pⁿ. Figure 2(b) depicts that all points in the level set function move in a horizontal direction when solving H–J equation. This condition can be used to explain an important drawback of CLSM that it cannot create new holes inside, i.e., lacking nucleation capability. Figure 3 illustrates several cases of updating of a level set surface, which means it can become steeper (Case 1) or flatter (Case 2) than before. Case 3 can never occur because the top part of the level set function cannot move downward to make a concave pit, as shown in Fig. 3. Case 4, in which the level set function is “pulled up”, also cannot happen due to the same reason. Figure 3 demonstrates a theoretical situation, but in practical numerical implementations, \mathrm{\Delta }t is not an infinitesimal value. The numerical errors may occasionally result in Cases 3 and 4 if the re-initialization scheme is not applied to recreate the signed distance level set function frequently.

Fig.2 Updating scheme for the level set function by solving Hamilton–Jacobi equation. (a) Two components for updating Φ; (b) all points move horizontally.

Full size|PPT slide

As mentioned above, CLSM is inherently lacking in the ability of nucleation capacity and is often criticized because it can induce the initial design-dependent problem, which means the final optimal design may heavily rely on the initial design. In the literature, researchers include numerous holes in the initial design to alleviate this problem. No widely recognized method exists for guiding people on how to establish the initial design. This drawback can be overcome by incorporating other nucleation schemes, such as topological derivatives [¹³], but CLSM still hardly handles topological changes naturally as density-based methods do.

Fig.3 Possible and impossible updated level set functions by solving Hamilton–Jacobi equation.

Full size|PPT slide

Zero level set method

When CLSM was introduced into the topology optimization field, a set of variational LSMs, such as a method based on dynamic implicit surface function [³⁶,³⁷], piecewise constant LSM [³⁸], LSM with a reaction–diffusion equation [³⁹], the topology representation function [⁴⁰], parameterized LSM [⁴¹–⁴⁵], and other methods [⁴⁶], was developed to overcome the drawbacks of CLSM. Most of these methods can be considered “LSM”, but they do not follow the solution procedure of CLSM. We suggest that these variational LSMs, which also use the zero level set to represent the boundary of design but do not exactly follow the set of conventional numerical manipulations to update the level set function, can be called “ZLSMs” to distinguish them from CLSM. These ZLSMs are all laudable attempts to explore the practical implementations of LSMs and provide valuable experiences and references for the endless further explorations.

The key point of these ZLSMs is that they also use a clear boundary to represent a design during the optimization process. Nonetheless, the numerical operation for level set updating is different. A way to realize ZLSMs is to revise Eq. (2) into Eq. (3) by directly removing the spatial difference term \left|\nabla \mathrm{\Phi }\right|:

(3){\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}t}}-V_{\mathrm{n}}=\mathrm{0.}

Thus, the updating of the level set function \mathrm{\Phi } only depends on the boundary normal velocity V_{\mathrm{n}}. The original H–J equation becomes an ordinary differential equation. The updating scheme of the revised Eq. (3) can be illustrated as Fig. 4.

Fig.4 Updating scheme of a zero level set method. (a) Updating scheme of Φ for Eq. (3); (b) all points move vertically.

Full size|PPT slide

In accordance with Fig. 4, V_{\mathrm{n}} indicates the vertical velocity of each point on the level set surface. The physical meaning of V_{\mathrm{n}} has changed, and the notation V_{\mathrm{n}} becomes unsuitable for it. The meaning of V_{\mathrm{n}} in this model is similar to the sensitivity in density-based methods. The calculation of V_{\mathrm{n}} can be borrowed from SIMP or ESO/BESO models. If \left|\nabla \mathrm{\Phi }\right|\mathrm{=1} is held with re-initialization as CLSM always executes, the effects of V_{\mathrm{n}} in those models become the same. We keep using notation V_{\mathrm{n}} in Eq. (3) for convenience to compare Eqs. (2) and (3). Although the inherent updating logic of Eqs. (2) and (3) relatively differ in accordance with Figs. 2 and 4, their convergence conditions are equivalent, i.e., V_{\mathrm{n}} equals zero due to the nonzero property of \left|\nabla \mathrm{\Phi }\right| along the boundary. Therefore, we can use the same V_{\mathrm{n}} to update the level set function in CLSM and ZLSM.

In our opinion, ZLSM may be more suitable for topology optimization than CLSM. The reasons that support our opinion are as follows. First, ZLSM naturally has nucleation capability, which can greatly alleviate the initial design-dependent problems, compared with CLSM. Second, ZLSM is more adaptable to miscellaneous meshes than CLSM. ZLSM uses an ordinary differential equation to update the level set function and has less complexity than the partial differential equation-driven CLSM, which is more convenient to be solved in a structured grid. Thus, ZLSM can be easily implemented in practical engineering problems in complex design domains and with unstructured meshes. A reasonable way in a commercial software system can be easily realized by sharing the same set of mesh with finite element analysis (FEA). The sensitivity analysis of ZLSM for certain problems can also be directly borrowed from density-based methods, and this point can make it easier to be implemented and accepted by people, such as for topology optimization on curved shell structures and sophisticated mathematical tools are needed for CLSM [⁴⁷]. The characteristics of the “customized” LSM in topology optimization can be summarized as it should be well connected to FEM and can be implemented as easily as the SIMP or ESO/BESO methods. Its advantages can be further developed.

Although ZLSM overcomes certain drawbacks of CLSM, the boundary evolution-based strategy in the level set-based model remains less natural than density-based methods in dealing with topological changes. In density-based methods, the topological optimization problem is transformed into a size optimization problem, and the topological change becomes a continuous process [³,⁴]. In the level set-based method, the topological change is a discrete process. The objective and constraint functions in the feasible space are not as continuous as in the density-based methods, which involve more difficulties in handling the topology optimization. This condition may also be considered a limitation of LSM; consequently, it is difficult to be commercialized. At this point, the density-based methods have predominant advantages in topology optimization. This issue has inspired us to introduce a density interpolation scheme to provide similar continuity in the level set-based method. In this study, we propose a new method that can improve the continuity of objective and constraint functions by using the high-dimensional information of the level set function. This method provides a simple way to combine the advantages of density-based and level set-based methods for realizing a natural topology evolution as a density-based method and a clear boundary representation solution.

Level set band method

In this section, a level set band method, which can be considered a combination of LSM with a density-based method to utilize the advantages of both methods, is proposed. This method follows, but is not restricted by, the parameterized LSM [⁴³] by incorporating a new parameter, i.e., level set band {\mathrm{\Phi }}_{\mathrm{b}}, which is the distance between a user-defined upper bound {\mathrm{\Phi }}_{\mathrm{u}} and a lower bound {\mathrm{\Phi }}_{\mathrm{l}}.

(4){\mathrm{\Phi }}_{\mathrm{b}}={\mathrm{\Phi }}_{\mathrm{u}}-{\mathrm{\Phi }}_{\mathrm{l}}.

The upper bound {\mathrm{\Phi }}_{\mathrm{u}} and lower bound {\mathrm{\Phi }}_{\mathrm{l}} do not mean the maximum and minimum values of the level set function, respectively, but indicate a range between which the densities of the elements should be interpolated with the values of the level set function. In this method, the density of each element in the structure and sensitivity analysis depend on its nodal values of the level set function. In the numerical implementation, the level set function value {\mathrm{\Phi }}_i in the middle of the ith element by interpolation can be adopted. Thus, the density of the ith element {\rho }_i can be calculated as

(5){\rho }_i=\mathrm{H}\left({\mathrm{\Phi }}_i\right),

where \mathrm{H}\in \{f:\mathrm{R}\to {\mathrm{R}}^{+}\} is the Heaviside function; here, it is numerically approximated as [¹⁴]

(6)

where ε is a small value indicates the lower bound density of void material and we can define \mathrm{\Delta }={\mathrm{\Phi }}_{\mathrm{u}}=-{\mathrm{\Phi }}_{\mathrm{l}} and {\mathrm{\Phi }}_{\mathrm{b}}=\mathrm{2}\mathrm{\Delta }.

Figure 5 presents the schematic of the proposed method. Figure 5(a) shows the interpolated density distribution for FEA with the level set function when {\mathrm{\Phi }}_{\mathrm{b}} is adequately large, and Fig. 5(b) illustrates the density distribution when {\mathrm{\Phi }}_{\mathrm{b}} is smaller. Thus, the design is determined by a narrow band boundary. If {\mathrm{\Phi }}_{\mathrm{b}} decreases to zero, the design is then determined by the zero level set, where \mathrm{\Phi }=\mathrm{0}. The optimization model becomes a ZLSM.

In the proposed method, the strategy is to gradually reduce the width of the level set band {\mathrm{\Phi }}_{\mathrm{b}} during the optimization process, as shown in Fig. 5. When \mathrm{\Delta } is sufficiently large, as shown in Fig. 5(a), all of the level set function values fall into the band {\mathrm{\Phi }}_{\mathrm{b}}. The density of each element is thus a projection of the related level set function value, similar to a definition of a density-based method, such as SIMP, where the design variables are defined as the density \rho and the upper and lower bounds are fixed at 1 and 0 (usually, a small value \epsilon instead), respectively. The strategy in SIMP is that the density of each element is forced to approach the upper or lower bounds (1 and 0) with a penalization scheme. Then, a black-and-white optimal design can be obtained. Unlike SIMP, the proposed method has variable upper and lower bounds in mapping the level set function. The distance between the upper and lower bounds is defined as a density band {\mathrm{\Phi }}_{\mathrm{b}}=\mathrm{2}\mathrm{\Delta }, which is represented in Eq. (4), and it is reduced gradually during the optimization process. When \mathrm{\Delta } becomes a considerably small value or zero, as shown in Fig. 5(b), the design is defined by the zero level set, and only the densities of the elements in an excessively small region around the boundary need to be calculated with the projection function Eq. (5). This method becomes an LSM. The density-based method in the initial stage finally converges to an LSM by assigning the parameter \mathrm{\Delta }={\mathrm{\Phi }}_{\mathrm{b}}/\mathrm{2} to a large value to start the optimization with a density-based method and then gradually decreasing the value of \mathrm{\Delta } during the optimization iteration process. We can combine the density-based method and LSM by involving only one parameter \mathrm{\Delta } to utilize the advantages of both methods. The density-based method presents flexibility in topology change and is minimally dependent on the initial design; LSM has clearly defined boundaries.

Fig.5 Density interpolation scheme of the level set band method by gradually reducing {\mathrm{\Phi }}_{\mathrm{b}}. (a) Case with a large level set band between the upper and lower bounds in the beginning stage; (b) case with a small level set band after convergence.

Full size|PPT slide

A general way to implement the level set-based method can be considered the case with a small band between the upper and lower bounds, as shown in Fig. 5(b). LSM is usually applied to a fixed Euler grid. The mesh and boundary are difficult to be conformed, and the elements around the boundary are usually cut through by the zero level set. An accurate way to calculate the contribution of the stiffness of “half” elements is using the extended FEM [⁴⁸–⁵⁰]; another simpler but with lower accuracy way is to use an approximate density to represent the stiffness contribution of the element around the boundary [⁴³]. The latter can be considered a means to interpolate the density of elements falling into the band or cut by the boundaries, as illustrated in Fig. 5(b), where the level set band is defined as a very small value. The proposed method can be considered an extension of this scheme to a wide space range around the boundary and to the entire optimization process over time. This method can be implemented in a simple manner by implementing a slight modification on the 88-line MATLAB code of parameterized LSM with radial basis functions [⁴³]. However, this method is not limited to the parametric way but can also be used to the discrete way in the level set framework. The proposed level set band method not only can be combined with ZLSM, as demonstrated in this paper, but also can be implemented in the CLSM framework. Compared with ZLSM, the proposed method needs further efforts in handling complex design problems.

Fig.6 Illustration of the level set band method with an intermediate optimal example (\mathrm{\Delta }={\mathrm{\Phi }}_{\mathrm{b}}/\mathrm{2}=\mathrm{3.5}).

Full size|PPT slide

The level set band decreases from a large value to a small value. The proposed model reveals a gradual change from a density-based method to a level set-based method. An example of an intermediate solution is provided to illustrate the level set band method in Fig. 6, in which the level set band {\mathrm{\Phi }}_{\mathrm{b}}=\mathrm{2}\mathrm{\Delta }=\mathrm{7} (the length of the finite element is 1), and the density distribution is obtained on the basis of the level set surface. The level sets on the upper and lower bounds are also plotted. They approach the same design as the zero level set after the level set band {\mathrm{\Phi }}_{\mathrm{b}} decreases to zero. From the viewpoint of the level set-based model, the proposed method can be considered an extended LSM by replacing the zero level set with a variable level set band. From another viewpoint, the level set band method can be considered a variation of a density-based method by replacing the fixed band of the density distribution between 0 and 1 with an alterable value. In this model, penalization and filtering mechanisms, as in SIMP, need not be introduced; related numerical issues [⁵¹,⁵²] no longer need to be processed. Figure 7 provides a schematic explanation of this point and shows how the proposed method can be considered a combination of the density-based method and LSM.

Fig.7 Level set band method can be considered a variation of the density-based method and the LSM by replacing the fixed bands {\mathrm{\Phi }}_{\mathrm{b}} with an alterable one.

Full size|PPT slide

Similar concepts by changing the level set isosurface were also evaluated in previous studies. In the moving isosurface threshold method, which was proposed in Ref. [¹⁶], a variational isosurface threshold for response surface, which has a similar definition to that of the parameter \mathrm{\Delta }={\mathrm{\Phi }}_{\mathrm{b}}/\mathrm{2} in this study for the level set surface, is applied to adjust the boundary design and determined by the Karush–Kuhn–Tucker condition. The level set-based method [⁵³] uses different layers of a unique characteristic level set function to represent different designs for optimizing connectable graded microstructures. This concept can also be used in minimum distance control [⁵⁴]. In the SIMP-based model, different density projections can be obtained to conduct robust topology optimization by implementing different threshold values [⁵⁵]; this process is similar to applying different levels in the density distribution function [⁵⁶]. The BESO-based model [⁵⁷] also involves a level set isosurface, which is iteratively determined by the upper and lower bounds. The two bounds approach almost the same value after convergence. The smooth boundaries are always clearly determined to separate the design domain into black and white parts during the optimization process, and the intermediate density elements only occur along the boundaries on the elements that are cut by the zero level set. The numerical examples indicate that the method illustrates similar properties to those of ZLSM because it adopts the zero level set to represent the design, although it comes from the BESO method. All of these successful implementations of the level set concept demonstrate its tremendous potential; it deserves further development in a future study.

Optimization scheme

The optimization model of the level set-based method can be defined as

(7)\begin{array}{l}{\min }_{\mathrm{\Phi }}\mathrm{}J(\mathit{u}\mathrm{,}\mathrm{\Phi }\mathrm{)=}{\displaystyle {\int }_{\mathrm{D}}f\left(\mathit{u}\right)H\left(\mathrm{\Phi }\right)}\mathrm{d}\Omega \hfill \\ \mathrm{s}\mathrm{.t}\mathrm{.}\{\begin{array}{l}a(\mathit{u},\mathrm{\hspace{0.17em}}\mathit{v},\mathrm{\hspace{0.17em}}\mathrm{\Phi })=l(\mathit{v},\mathrm{\hspace{0.17em}}\mathrm{\Phi }),\mathrm{}\forall \mathit{v}\in \mathit{U},\hfill \\ G\left(\mathrm{\Phi }\right)={\displaystyle \frac{{\displaystyle {\int }_{D}H\left(\mathrm{\Phi }\right)}\mathrm{d}\Omega }{{\displaystyle {\int }_D\mathrm{d}\Omega }}}-V_{\max }\le \mathrm{0},\hfill \\ \mathit{u}={\mathit{u}}_{\mathrm{0}}\mathrm{in}{\mathrm{\Gamma }}_{\mathit{u}},\hfill \\ \mathit{C}:\mathit{\epsilon }\left(\mathit{u}\right)·\mathit{n}=\mathit{\tau }\mathrm{in}{\mathrm{\Gamma }}_{\mathit{\tau }},\hfill \end{array}\hfill \end{array}

where J is the objective function for a specific physical type described by f, u is the displacement field, \mathit{\epsilon } is the linearized strain tensor, \mathit{C} is the elasticity tensor, v is the adjoint displacement in the space \mathit{U} of the kinematically admissible displacement fields, G\left(\mathrm{\Phi }\right) is the volume constraint to limit material usage, V_{\max } is the maximum allowable volume fraction of the design domain, {\mathit{u}}_{\mathrm{0}}, \mathit{n}, and \mathit{\tau } are the given displacement, the boundary unit normal vector and traction, respectively, and {\mathrm{\Gamma }}_u and {\mathrm{\Gamma }}_{\tau } are the Dirichlet and Neumann boundaries, respectively. The energy bilinear form a\left(\mathit{u},\mathrm{\hspace{0.17em}}\mathit{v},\mathrm{\hspace{0.17em}}\mathrm{\Phi }\right) and the load linear form l\left(\mathit{v},\mathrm{\hspace{0.17em}}\mathrm{\Phi }\right) are defined as

(8)a(\mathit{u},\mathrm{\hspace{0.17em}}\mathit{v},\mathrm{\hspace{0.17em}}\mathrm{\Phi })={\displaystyle {\int }_D\left(\mathit{\epsilon }\left(\mathit{u}\right):\mathit{C}:\mathit{\epsilon }\left(\mathit{v}\right)\right)\mathrm{H}\left(\mathrm{\Phi }\right)\mathrm{d}\mathrm{\Omega }},

(9)\mathit{l}(\mathit{v},\hspace{0.17em}\mathrm{\Phi })={\displaystyle {\int }_{{\mathrm{\Gamma }}_{\mathit{\tau }}}\mathit{\tau }·\mathit{v}\mathrm{d}\mathrm{\Gamma }+}{\displaystyle {\int }_D\mathit{b}·\mathit{v}}\mathrm{H}\left(\mathrm{\Phi }\right)\mathrm{d}\mathrm{\Omega },

where \mathit{b} represents the body force.

In the framework of CLSM, the sensitivity analysis based on shape derivative [¹³,¹⁴] can be used to derive the normal velocity V_{\mathrm{n}} along the moving boundary in the steepest descent direction.

V_{\mathrm{n}}=-f\left(\mathit{u}\right)-\mathit{\epsilon }\left(\mathit{u}\right):\mathit{C}:\mathit{\epsilon }\left(\mathit{v}\right)

(10)+\mathit{b}·\mathit{v}+\left(\nabla \left(\mathit{\tau }·\mathit{v}\right)·\mathit{n}\right)+\kappa \left(\mathit{\tau }·\mathit{v}\right)-\lambda ,

where \kappa =\nabla ·\mathit{n} is the curvature around the boundary, λ is the Lagrange multiplier to control the volume constraint, and it can be calculated with the augmented Lagrange multiplier method [⁴³,⁵⁸] or bisectional method [²¹]. For a compliance minimization problem f\left(\mathit{u}\right)=\mathit{\epsilon }\left(\mathit{u}\right):\mathit{C}:\mathit{\epsilon }\left(\mathit{u}\right) without body force and boundary traction, the velocity can be simplified as

(11)V_{\mathrm{n}}=\mathit{\epsilon }\left(\mathit{u}\right):\mathit{C}:\mathit{\epsilon }\left(\mathit{u}\right)-\lambda .

With the velocity V_{\mathrm{n}} obtained in Eq. (11), the CLSM of Eq. (2) or ZLSM of Eq. (3) can be used to update the level set function until convergence to realize the topology optimization. This velocity is meaningful only along the boundary in CLSM because only the variation on the boundary can affect the objective function. In the proposed level set band method (Fig. 8), the evolution of the level set function of the solid part does not influence the objective function and volume constraint given that the changed part does not fall into the level set band, as shown in Fig. 8(b). The sensitivity with respect to the level set function can be considered zero, but updating the level set function at that interior area has the potential to change the objective function and the topology of the design, as shown in Fig. 8(c). This diagram also illustrates the scheme that the existing level set band makes the topology change a continuous procedure. In CLSM or ZLSM, the topology change is always a discontinuous procedure and may cause numerical issues, such as oscillation or local optimal.

Fig.8 Two cases of variation in the level set function without and with changing the objective function. (a) Original level set function; (b) the objective function is unchanged (Case 1); (c) the objective function is changed (Case 2).

Full size|PPT slide

In this model, the velocity V_{\mathrm{n}} is extended to the entire design domain by directly calculating the velocity with Eq. (11) over the design domain. This process can be considered a “natural velocity extension” approach applied in most level set-based topology optimization models [¹⁴,⁴¹]. Thus, the level set function can be updated in the entire design domain, and the nucleation capability is easily realized. This approach is unlike the conventional way to extend the velocity in CLSM to keep the level set function a signed distance function as applied in the computer graphics field [¹⁰,³⁴].

On the basis of Eq. (3), the level set function can be updated based on ZLSM with the first-order difference scheme as

(12){\mathrm{\Phi }}^{i+\mathrm{1}}={\mathrm{\Phi }}^i+V_{\mathrm{n}}^i·\mathrm{\Delta }t,

where the superscripts i and i+1 indicate the iteration steps, and \mathrm{\Delta }t is the time step size. In this study, the evolution of the level set function is realized by updating the coefficients in the parameterized LSM [⁴³]. Readers can refer to Ref. [⁴³] and the provided code for a detailed implementation approach. This scheme can also be implemented on the basis of CLSM by updating the level set function with Eq. (2) to obtain an accurate boundary evolution solution. Nevertheless, the related numerical manipulations of CLSM should be adopted, and the nucleation capability will be missing.

Numerical examples

In this section, several examples are analyzed to illustrate the effects of the proposed method. Please note the numerical examples in this paper are dimensionless. The basic parameters follow the 88-line MATLAB code [⁴³], including Young’s elasticity modulus E=\mathrm{1} for solid material, E={\mathrm{10}}^{-\mathrm{9}} for void material, and Poisson’s ratio \nu =\mathrm{0.3}. The used radial basis functions are multiquadratic spline with c={\mathrm{10}}^{-\mathrm{4}}. Only the mean compliance minimization problem is studied, but this method can be easily applied to other problems, such as compliance mechanism design and material design problems. The convergence condition is set to be the same as in the 88-line MATLAB code that in all of the last nine steps, the mean compliance M satisfies the following criterion:

(13)

where the subscript of M means the iterative step number.

Cantilever beam

The first example is shown in Fig. 9. A short cantilever beam is given, in which the left side is fixed. Concentrated force F=\mathrm{1} is applied vertically downward at the middle point at the right side. The size of the design domain is 80 × 40, and 80 × 40 four-node bilinear square elements are used to perform FEA. The total volume fraction is set as 50%.

Fig.9 Design domain and boundary conditions of a cantilever beam.

Full size|PPT slide

Figure 10 shows the iterative designs of the structure in Steps 1, 10, 20, 30, 50, and 109 (converged). Figure 10(a) shows the zero level set of the designs, Fig. 10(b) depicts the density distributions of the designs, Fig. 10(c) illustrates the level sets of the upper and lower bounds (red is the upper bound, and blue is the lower bound), and Fig. 10(d) indicates the level set functions and the upper and lower bound planes during the optimization. The initial value of \mathrm{\Delta } is 5 and decreases by 0.1 at each step to the minimum value of 0.5 after 45 steps.

The initial value of the level set function \mathrm{\Phi } is given as -\mathrm{3}\le \mathrm{\Phi }\le \mathrm{3}, and the initial design is filled with gray elements, as shown in Fig. 10(b). In the first figure of Fig. 10(c), the initial level set function does not touch the upper and lower bounds (\pm \mathrm{5}), and thus, no red or white area exists. After several iterations, the values of the level set function and Δ change, as shown in the second figure in Fig. 10(d). The red part in Fig. 10(c) indicates the area where the density is 1, and the white part indicates the area where the density is ε. Figure 10(a) indicates the zero level set, but it is not the real design because the real finite element model for analysis is given by Fig. 10(b). This problem converges after 109 iterations and then the final zero level set becomes consistent with the finite element model, as depicted in Figs. 10(a) and 10(b). The level sets on the upper and lower bound planes almost perfectly coincide because the distance between the upper and the lower bounds is very small ({\mathrm{\Phi }}_{\mathrm{b}}=\mathrm{2}\mathrm{\Delta }=\mathrm{1}).

The iteration process of this example can be considered a model of the density-based method; if the initial value of \mathrm{\Delta } is set to a small value (e.g., 0.5), the iteration process then becomes a zero level set model. The results are given in Fig. 11, and the iteration numbers are 1, 10, 20, 30, 50, and 131. Figure 11(b) depicts that the design has few gray elements around the boundary from the initial design to the final design. The zero level set is almost the same as the density distribution. Figure 11(c) also shows that the level sets on the upper and lower bound planes are almost coincident from the start to the end. Compared with the optimization process shown in Fig. 10, this case can be considered a zero level set model because the clear boundary is always given. Therefore, boundary-related problems, such as gravity or hydraulic pressures, can be easily handled.

Fig.10 Optimization iteration process of the 2D cantilever beam with decreasing Δ from 5 to 0.5 by 0.1 at each step (from top to bottom: Steps 1, 10, 20, 30, 50, and 109). (a) Zero level set; (b) density distribution; (c) level sets on the lower and upper bound planes; (d) level set function and the lower and upper bound planes.

Full size|PPT slide

Fig.11 Optimization iteration process of the 2D cantilever beam with fixed \mathrm{\Delta }=\mathrm{0.5} (from top to bottom: Steps 1, 10, 20, 30, 50, and 131). (a) Zero level set; (b) density distribution; (c) level sets on the lower and upper bound planes; (d) level set function and the lower and upper bound planes.

Full size|PPT slide

In the following part, the same problem is solved with different initial designs to illustrate the initial design-dependent issue of the proposed method. Figures 12(a) and 12(b) show the case with decreasing Δ from 5 to 0.5 by 0.1 at each step, and Figs. 12(c) and 12(d) show the case with fixed \mathrm{\Delta }\mathrm{=0}\mathrm{.5}. Figures 12(a) and 12(c) are the zero level sets during the iteration, and Figs. 12(b) and 12(d) are the corresponding density distributions. In the two cases, the initial density of each element is 0.5.

Fig.12 Optimization process of the cantilever beam problem with uniformly distributed initial density. (a) Zero level set: Δ decreases from 5 to 0.5 by 0.1 at each step, and the iteration steps are 1, 10, 20, 30, 50, and 149 (from top to bottom); (b) density distribution: Δ decreases from 5 to 0.5 by 0.1 at each step, and the iteration steps are 1, 10, 20, 30, 50, and 149 (from top to bottom); (c) zero level set: fixed \mathrm{\Delta }=\mathrm{0.5}, and the iteration numbers are 1, 10, 20, 30, 50, and 126 (from top to bottom); (d) density distribution: fixed \mathrm{\Delta }=\mathrm{0.5}, and the iteration numbers are 1, 10, 20, 30, 50, and 126 (from top to bottom).

Full size|PPT slide

Figure 12 demonstrates that the density distributions and topologies in the intermediate steps of the two cases are relatively different. In Figs. 12(a) and 12(b), the topology change is driven by the density variation in the earlier stage, in which Figs. 12(a) and 12(b) can be considered a density-based model. In Figs. 12(c) and 12 (d), the density distribution is almost black and white everywhere during the optimization process; this can be considered a zero level set model with nucleation. The only difference between the two methods is the value of Δ or {\mathrm{\Phi }}_{\mathrm{b}}. The comparison in this example clearly illustrates the most important characteristic of the proposed model, i.e., the density-based method and LSM can be connected with only one parameter Δ. If Δ decreases from a large value to a small value, this model becomes a density-based model; if Δ is set as a small value from the beginning to the end, this model becomes a zero level set model. The penalization and density filter are not needed to be specifically applied in this model.

Different Δ updating schemes for a simply supported beam problem

In the following part, a simply supported beam, shown in Fig. 13, is evaluated with the same method to illustrate the convergence property of the proposed method. The beam size is 160 × 40 and is discretized to 160 × 40 Q4 elements. On the basis of the symmetry of the problem, only the right half part of the design is optimized.

Fig.13 Design domain and boundary conditions of a simply supported beam.

Full size|PPT slide

The simply supported beam problem generally varies in terms of the topology of the final designs. In this part, different initial designs and three Δ updating Schemes (a), (b), and (c) are imposed to provide a clearer understanding of the proposed model as described in Figs. 14(a), 14(b), and 14(c), respectively. Here, three Δ updating schemes are implemented. Scheme (a) uses a fixed Δ value of 0.01 in the overall optimization process; Schemes (b) and (c) use decreasing Δ values with different reductions. {\mathrm{\Delta }}_{\mathrm{0}} denotes the initial Δ value, \mathrm{d}\mathrm{\Delta } denotes the reduction value of Δ in each iteration, and \min \mathrm{\Delta } denotes the minimum value of Δ. N indicates the total iteration number when the optimization convergences, and MC means the mean compliance of the design.

Fig.14 Optimal results of the simply supported beam with different initial designs and three schemes. (a) \mathrm{\Delta }\mathrm{=}{\mathrm{\Delta }}_{\mathrm{0}}\mathrm{=}\min \mathrm{\Delta }=\mathrm{0.01}; (b) {\mathrm{\Delta }}_{\mathrm{0}}=\mathrm{5}, \mathrm{d}\mathrm{\Delta }=-\mathrm{0.1}, \min \mathrm{\Delta }=\mathrm{0.01}; (c) {\mathrm{\Delta }}_{\mathrm{0}}=\mathrm{5}, \mathrm{d}\mathrm{\Delta }=-\mathrm{0.05}, \min \mathrm{\Delta }=\mathrm{0.01}.

Full size|PPT slide

As shown in Fig. 14, in each problem, the upper figure is the initial design, and the lower figure is the final design after convergence. In each updating scheme, three randomly generated initial designs are used to examine the algorithm, and the last one is a fixed initial design with three holes inside. The iteration numbers and objective functions are provided in Fig. 14. These data are plotted in Fig. 15, in which the distribution of the final results can be clearly studied.

From the previous study, Scheme (a) can be considered a ZLSM. Based on the results shown in Figs. 14 and 15, it can be found with the random initial design, the results show that this scheme needs relatively more iteration steps and obtains high objective function values. When given an appropriate initial design, the scheme obtains the best result. This observation shows that ZLSM generally needs a proper initial guess; otherwise, it may not achieve a good result. This condition is considered the initial design-dependent problem. The nucleation capability can greatly alleviate this problem, but the initial guess is still important for level set-based models.

Schemes (b) and (c) involve the density distribution stage that makes them close to a density-based optimization model. Scheme (b) uses a larger reduction of Δ in each iteration and shows faster convergent speed compared with those of Scheme (c). A gradual decrease in Δ can generate a reasonable final design, as demonstrated in Fig. 15; all of the final designs using Scheme (c) have low final mean compliances. Another observation is that the initial design for the density-involved schemes is also important because it accelerates the convergence greatly in Scheme (b). In Scheme (c), the results are adequately good, and the improvement is thus not obvious.

This numerical study can be concluded with the following suggestion: If an appropriate initial design can be easily obtained, Schemes (a) and (b) should be chosen; otherwise, a small value for decreasing the level set band as in Scheme (c) should be selected to obtain a reasonable design by sacrificing some efficiency.

Fig.15 Comparison of the objective functions and the total iteration steps of the simply supported beam problem with different initial designs and Δ values.

Full size|PPT slide

3D cantilever beam

This method can also be applied to 3D models without difficulties. Figure 16 illustrates the optimization iteration process of a 3D cantilever beam optimization problem. The left side of the beam is fixed, and a downward force is applied at the middle point of the bottom line of the right side. The structure is discretized with 60 × 30 × 10 elements. The volume fraction is set at 20%. Δ is set at 5 at the start and decreased by 0.1 each step to 0.5. Figure 16(a) is the zero level set of the design at Steps 1, 15, 30, and 130. The corresponding density distributions are shown in Fig. 16(b) by the way provided by Ref. [⁵⁹]. The initial design has a uniformly distributed density of 0.5. The iteration process can be considered a density-based model. In the end, only a small number of elements are gray around the boundary. The iteration process illustrates that the topology change can be easily realized, and a clear boundary can be obtained without implementing any penalization and filter schemes.

Fig.16 Implementation of the proposed method in the optimization of a 3D cantilever beam. (a) Zero level set of the 3D cantilever beam problem (left to right: Steps 1, 15, 30, and 130); (b) density distribution of the 3D cantilever beam problem (left to right: Steps 1, 15, 30, and 130).

Full size|PPT slide

Conclusions

In this paper, a comparison between CLSM and ZLSM for solving structural topology optimization problems is schematically discussed. ZLSM, which can be easily applied and solved with minimal numerical issues, is suggested in practical implementations. We propose the level set band method to learn from density-based methods with a remarkable flexibility change in topology and improve the level set-based method by introducing a level set band {\mathrm{\Phi }}_{\mathrm{b}}, which utilizes the high-dimensional information of the level set function to improve the continuity of the objective and constraint functions in handling topology changes. The density and the level set-based methods can be seamlessly combined by changing the value of {\mathrm{\Phi }}_{\mathrm{b}} gradually. The proposed model with a large value of the level set band illustrates the property of the density-based model and is highly flexible in handling topology change. When the size of the level set band is decreased, this method becomes similar to an LSM, which can provide clear boundaries or a black-and-white design, without involving penalization. Thus, the parameter, level set band, can be used to connect the density-based method and LSM. In another aspect, this algorithm shows that the two methods have no essential difference. Numerical examples with random initial designs are used to evaluate the convergence property of the proposed method. If the initial design is unclear, the level set band can be slowly decreased to obtain a highly reasonable design with lesser efficiency; otherwise, the size of the level set band can be decreased rapidly to accelerate the convergence process. 2D and 3D examples are solved to illustrate the effectiveness of the proposed method.

Acknowledgments

The authors gratefully acknowledge the financial support provided by the National Natural Science Foundation of China (Grant No. 11372004), the State Key Laboratory of Subtropical Building Science (Grant No. 2016 KB13), and the State Key Laboratory of Structural Analysis for Industrial Equipment (Grant No. GZ18109).

Open Access

This article is licensed under Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution, and reproduction in any medium or format, as long as appropriate credit is given to the original author(s) and the source, a link is provided to the Creative Commons license, and any changes made are indicated.ƒImages or other third-party materials in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license 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 license, visit http://creativecommons.org/licenses/by/4.0/.

References

1 Cheng K T, Olhoff N. An investigation concerning optimal design of solid elastic plates. International Journal of Solids and Structures, 1981, 17(3): 305–323

2 Bendsøe M P, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 197–224

3 Bendsøe M P. Optimal shape design as a material distribution problem. Structural Optimization, 1989, 1(4): 193–202

4 Rozvany G I, Zhou M, Birker T. Generalized shape optimization without homogenization. Structural Optimization, 1992, 4(3–4): 250–252

5 Bendsøe M P, Sigmund O. Material interpolation schemes in topology optimization. Archive of Applied Mechanics, 1999, 69(9–10): 635–654

6 Stolpe M, Svanberg K. An alternative interpolation scheme for minimum compliance topology optimization. Structural and Multidisciplinary Optimization, 2001, 22(2): 116–124

7 Xie Y M, Steven G P. A simple evolutionary procedure for structural optimization. Computers & Structures, 1993, 49(5): 885–896

8 Querin O M, Steven G P, Xie Y M. Evolutionary structural optimisation (ESO) using a bidirectional algorithm. Engineering Computations, 1998, 15(8): 1031–1048

9 Huang X, Xie Y M. Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method. Finite Elements in Analysis and Design, 2007, 43(14): 1039–1049

10 Osher S, Sethian J A. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics, 1988, 79(1): 12–49

11 Sethian J A, Wiegmann A. Structural boundary design via level set and immersed interface methods. Journal of Computational Physics, 2000, 163(2): 489–528

12 Osher S J, Santosa F. Level set methods for optimization problems involving geometry and constraints I. Frequencies of a two-density inhomogeneous drum. Journal of Computational Physics, 2001, 171(1): 272–288

13 Allaire G, Jouve F, Toader A M. Structural optimization using sensitivity analysis and a level-set method. Journal of Computational Physics, 2004, 194(1): 363–393

14 Wang M Y, Wang X M, Guo D M. A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1–2): 227–246

15 Sui Y K, Ye H L. Continuum Topology Optimization Methods ICM. Beijing: Science Press, 2013 (in Chinese)

16 Tong L Y, Lin J Z. Structural topology optimization with implicit design variable—Optimality and algorithm. Finite Elements in Analysis and Design, 2011, 47(8): 922–932

17 Wei P, Ma H T, Wang M Y. The stiffness spreading method for layout optimization of truss structures. Structural and Multidisciplinary Optimization, 2014, 49(4): 667–682

18 Cao M J, Ma H T, Wei P. A modified stiffness spreading method for layout optimization of truss structures. Acta Mechanica Sinica, 2018, 34(6): 1072–1083

19 Zhang W S, Yuan J, Zhang J, A new topology optimization approach based on moving morphable components (MMC) and the ersatz material model. Structural and Multidisciplinary Optimization, 2016, 53(6): 1243–1260

20 Zhang W S, Li D D, Kang P, Explicit topology optimization using IGA-based moving morphable void (MMV) approach. Computer Methods in Applied Mechanics and Engineering, 2020, 360: 112685

21 Sigmund O. A 99 line topology optimization code written in Matlab. Structural and Multidisciplinary Optimization, 2001, 21(2): 120–127

22 Huang X, Xie Y M. Evolutionary Topology Optimization of Continuum Structures: Methods and Applications. New York: Wiley, 2010

23 Osher S, Fedkiw R. Level Set Methods and Dynamic Implicit Surfaces. New York: Springer, 2002

24 Sethian J A. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science. London: Cambridge University Press, 1999

25 Xia Q, Shi T, Wang M Y, Simultaneous optimization of cast part and parting direction using level set method. Structural and Multidisciplinary Optimization, 2011, 44(6): 751–759

26 van Dijk N P, Maute K, Langelaar M, van Keulen F. Level-set methods for structural topology optimization: A review. Structural and Multidisciplinary Optimization, 2013, 48(3): 437–472

27 Xia Q, Shi T. Optimization of structures with thin-layer functional device on its surface through a level set based multiple-type boundary method. Computer Methods in Applied Mechanics and Engineering, 2016, 311: 56–70

28 Jiang L, Guo Y, Chen S K, Concurrent optimization of structural topology and infill properties with a CBF-based level set method. Frontiers of Mechanical Engineering, 2019, 14(2): 171–189

29 Xia Q, Shi T, Xia L. Stable hole nucleation in level set based topology optimization by using the material removal scheme of BESO. Computer Methods in Applied Mechanics and Engineering, 2019, 343: 438–452

30 Barrera J L, Geiss M J, Maute K. Hole seeding in level set topology optimization via density fields. Structural and Multidisciplinary Optimization, 2020, 61: 1319–1343

31 Shu L, Wang M Y, Fang Z D, Level set based structural topology optimization for minimizing frequency response. Journal of Sound and Vibration, 2011, 330(24): 5820–5834

32 Shu C W, Osher S. Efficient implementation of essentially nonoscillatory shock-capturing schemes. II. Journal of Computational Physics, 1989, 83(1): 32–78

33 Jiang G S, Peng D P. Weighted ENO schemes for Hamilton–Jacobi equations. SIAM Journal on Scientific Computing, 2000, 21(6): 2126–2143

34 Peng D P, Merriman B, Osher S, A PDE-based fast local level set method. Journal of Computational Physics, 1999, 155(2): 410–438

35 Adalsteinsson D, Sethian J A. The fast construction of extension velocities in level set methods. Journal of Computational Physics, 1999, 148(1): 2–22

36 Guo X, Zhao K, Wang M Y. A new approach for simultaneous shape and topology optimization based on dynamic implicit surface function. Control and Cybernetics, 2005, 34(1): 255–282

37 Belytschko T, Xiao S P, Parimi C. Topology optimization with implicit functions and regularization. International Journal for Numerical Methods in Engineering, 2003, 57(8): 1177–1196

38 Wei P, Wang M Y. Piecewise constant level set method for structural topology optimization. International Journal for Numerical Methods in Engineering, 2009, 78(4): 379–402

39 Yamada T, Izui K, Nishiwaki S, A topology optimization method based on the level set method incorporating a fictitious interface energy. Computer Methods in Applied Mechanics and Engineering, 2010, 199(45–48): 2876–2891

40 de Ruiter M J, van Keulen F. Topology optimization using a topology description function. Structural and Multidisciplinary Optimization, 2004, 26(6): 406–416

41 Wang S Y, Wang M Y. Radial basis functions and level set method for structural topology optimization. International Journal for Numerical Methods in Engineering, 2006, 65(12): 2060–2090

42 Luo Z, Wang M Y, Wang S Y, A level set-based parameterization method for structural shape and topology optimization. International Journal for Numerical Methods in Engineering, 2008, 76(1): 1–26

43 Wei P, Li Z Y, Li X P, An 88-line MATLAB code for the parameterized level set method based topology optimization using radial basis functions. Structural and Multidisciplinary Optimization, 2018, 58(2): 831–849

44 Wang Y J, Benson D J. Geometrically constrained isogeometric parameterized level-set based topology optimization via trimmed elements. Frontiers of Mechanical Engineering, 2016, 11(4): 328–343

45 Wei P, Paulino G H. A parameterized level set method combined with polygonal finite elements in topology optimization. Structural and Multidisciplinary Optimization, 2020, 61: 1913–1928

46 Burger M, Hackl B, Ring W. Incorporating topological derivatives into level set methods. Journal of Computational Physics, 2004, 194(1): 344–362

47 Ye Q, Guo Y, Chen S K, Topology optimization of conformal structures on manifolds using extended level set methods (X-LSM) and conformal geometry theory. Computer Methods in Applied Mechanics and Engineering, 2019, 344: 164–185

48 Wei P, Wang M Y, Xing X H. A study on X-FEM in continuum structural optimization using a level set model. Computer Aided Design, 2010, 42(8): 708–719

49 Li L, Wang M Y, Wei P. XFEM schemes for level set based structural optimization. Frontiers of Mechanical Engineering, 2012, 7(4): 335–356

50 Geiss M J, Barrera J L, Boddeti N, A regularization scheme for explicit level-set XFEM topology optimization. Frontiers of Mechanical Engineering, 2019, 14(2): 153–170

51 Bruyneel M, Duysinx P. Note on topology optimization of continuum structures including self-weight. Structural and Multidisciplinary Optimization, 2005, 29(4): 245–256

52 Zhu J H, He F, Liu T, Structural topology optimization under harmonic base acceleration excitations. Structural and Multidisciplinary Optimization, 2018, 57(3): 1061–1078

53 Wang Y Q, Chen F F, Wang M Y. Concurrent design with connectable graded microstructures. Computer Methods in Applied Mechanics and Engineering, 2017, 317: 84–101

54 Kang Z, Wang Y G, Wang Y Q. Structural topology optimization with minimum distance control of multiphase embedded components by level set method. Computer Methods in Applied Mechanics and Engineering, 2016, 306: 299–318

55 Guest J K, Prevost J H, Belytschko T. Achieving minimum length scale in topology optimization using nodal design variables and projection functions. International Journal for Numerical Methods in Engineering, 2004, 61(2): 238–254

56 Wang F W, Lazarov B S, Sigmund O. On projection methods, convergence and robust formulations in topology optimization. Structural and Multidisciplinary Optimization, 2011, 43(6): 767–784

57 Da D C, Xia L, Li G Y, Evolutionary topology optimization of continuum structures with smooth boundary representation. Structural and Multidisciplinary Optimization, 2018, 57(6): 2143–2159

58 Rockafellar R T. The multiplier method of Hestenes and Powell applied to convex programming. Journal of Optimization Theory and Applications, 1973, 12(6): 555–562

59 Liu K, Tovar A. An efficient 3D topology optimization code written in Matlab. Structural and Multidisciplinary Optimization, 2014, 50(6): 1175–1196