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

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.


Introduction
The successful topology optimization of continuums inspired from the study of the optimal thickness distribution of elastic plates [1] has led to the successive development of several topology optimization models, such as the homogenization method [2], solid isotropic material/microstructure with penalization (SIMP) method [3][4][5], rational approximation of material property method [6], evolutionary structural optimization (ESO)/bidirectional ESO (BESO) [7][8][9], level set method (LSM) [10][11][12][13][14], independent-continuous mapping method [15], moving isosurface threshold method [16], stiffness spreading method [17,18], and moving morphable component/void method [19,20]. 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 [21]. 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 [7], 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. [8,9,22] 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 [10][11][12][13][14][23][24][25][26][27][28][29][30], which is borrowed from computer graphics and image-processing fields, was first implemented in topology optimization around 2000 [11,12]. 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 [13,14]. 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 [31]. 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 [10] for interface tracking. It has been successfully implemented in computer graphics and image segmentation fields before introducing into the topology optimization field [23,24]. The concept of a 2D design, which is represented with a level set function, is illustrated in Fig. 1.
The basic idea of LSM is to define a one dimension higher surface, i.e., the level set function Φ, 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: where x 2 R 2 or R 3 denotes a point in the design domain D & R 2 or R 3 , Ω and ∂Ω 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 [23,24], 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 [32,33], the reinitialization process to recreate a signed distance level set function [24,34], and the velocity extension approach [35] 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 [23,24]; 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 rΦ and temporal derivative ∂Φ=∂t: where V 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 [14]. The " -" in Eq.
(2) depends on the outside positive normal direction of V n , which is determined by the level set representation model given in Eq. (1). If Φðx, tÞ > 0 is defined as the outside part, then " -" should be changed to "þ". 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 Δt. _ Φ in Eq. (2) consists of two parts, which are ð∂Φ=∂tÞΔt and V n jrΦjΔt. The two parts are illustrated in Fig. 2(a). The equation _ Φ ¼ 0 is always satisfied in each step. Therefore, the two parts are cancelled each other out, and a point on the level set surface Φ of the nth time step p n moves horizontally to point p n +1 in the next time step. Here, V n denotes the normal velocity on point p n . 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, Δ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.
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 [13], but CLSM still hardly handles topological changes naturally as density-based methods do.

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 [36,37], piecewise constant LSM [38], LSM with a reactiondiffusion equation [39], the topology representation function [40], parameterized LSM [41][42][43][44][45], and other methods [46], 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 jrΦj:  In accordance with Fig. 4, V n indicates the vertical velocity of each point on the level set surface. The physical meaning of V n has changed, and the notation V n becomes unsuitable for it. The meaning of V n in this model is similar to the sensitivity in density-based methods. The calculation of V n can be borrowed from SIMP or ESO/BESO models. If jrΦj ¼ 1 is held with re-initialization as CLSM always executes, the effects of V n in those models become the same. We keep using notation V 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 n equals zero due to the nonzero property of jrΦj along the boundary. Therefore, we can use the same V 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 [47]. 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 densitybased 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   [3,4]. 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 densitybased 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 [43] by incorporating a new parameter, i.e., level set band Φ b , which is the distance between a user-defined upper bound Φ u and a lower bound Φ l .
The upper bound Φ u and lower bound Φ 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 Φ i in the middle of the ith element by interpolation can be adopted. Thus, the density of the ith element i can be calculated as where H 2 ff : R↕ ↓ R þ g is the Heaviside function; here, it is numerically approximated as [14] HðxÞ ¼ where ε is a small value indicates the lower bound density of void material and we can define  Figure 5(a) shows the interpolated density distribution for FEA with the level set function when Φ b is adequately large, and Fig. 5(b) illustrates the density distribution when Φ b is smaller. Thus, the design is determined by a narrow band boundary. If Φ b decreases to zero, the design is then determined by the zero level set, where Φ ¼ 0. The optimization model becomes a ZLSM.
In the proposed method, the strategy is to gradually reduce the width of the level set band Φ b during the optimization process, as shown in Fig. 5. When Δ is sufficiently large, as shown in Fig. 5(a), all of the level set function values fall into the band Φ 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 and the upper and lower bounds are fixed at 1 and 0 (usually, a small value ε 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 Φ b ¼ 2Δ, which is represented in Eq. (4), and it is reduced gradually during the optimization process. When Δ 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 Δ ¼ Φ b =2 to a large value to start the optimization with a density-based method and then gradually decreasing the value of Δ during the optimization iteration process. We can combine the density-based method and LSM by involving only one parameter Δ 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.
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 [48][49][50]; another simpler but with lower accuracy way is to use an approximate density to represent the stiffness contribution of the element around the boundary [43]. 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 [43]. 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.
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 Φ b ¼ 2Δ ¼ 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 Φ b  decreases to zero. From the viewpoint of the level setbased 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 densitybased 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 [51,52] 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 densitybased method and LSM. 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. [16], a variational isosurface threshold for response surface, which has a similar definition to that of the parameter Δ ¼ Φ b =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 [53] 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 [54]. In the SIMP-based model, different density projections can be obtained to conduct robust topology optimization by implementing different threshold values [55]; this process is similar to applying different levels in the density distribution function [56]. The BESO-based model [57] 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 where J is the objective function for a specific physical type described by f, u is the displacement field, ε is the linearized strain tensor, C is the elasticity tensor, v is the adjoint displacement in the space U of the kinematically admissible displacement fields, GðΦÞ is the volume constraint to limit material usage, V max is the maximum allowable volume fraction of the design domain, u 0 , n, and τ are the given displacement, the boundary unit normal vector and traction, respectively, and Γ u and Γ τ are the Dirichlet and Neumann boundaries, respectively. The energy bilinear form aðu, v, ΦÞ and the load linear form lðv, ΦÞ are defined as aðu, v, ΦÞ ¼ ! D εðuÞ : C : εðvÞ HðΦÞdΩ, where b represents the body force.
In the framework of CLSM, the sensitivity analysis based on shape derivative [13,14] can be used to derive the normal velocity V n along the moving boundary in the steepest descent direction. Fig. 7 Level set band method can be considered a variation of the density-based method and the LSM by replacing the fixed bands Φ b with an alterable one.
where κ ¼ r$n is the curvature around the boundary, l is the Lagrange multiplier to control the volume constraint, and it can be calculated with the augmented Lagrange multiplier method [43,58] or bisectional method [21]. For a compliance minimization problem f ðuÞ ¼ εðuÞ : C : εðuÞ without body force and boundary traction, the velocity can be simplified as V n ¼ εðuÞ : C : εðuÞ -l: With the velocity V 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.
In this model, the velocity V 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 [14,41]. 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 [10,34].
On the basis of Eq. (3), the level set function can be updated based on ZLSM with the first-order difference scheme as where the superscripts i and i + 1 indicate the iteration steps, and Δ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 [43]. Readers can refer to Ref. [43] 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 [43], including Young's elasticity modulus E ¼ 1 for solid material, E ¼ 10 -9 for void material, and Poisson's ratio ¼ 0:3. The used radial basis functions are multiquadratic spline with c ¼ 10 -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: 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 ¼ 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%. 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 Δ 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 Φ is given as -3£Φ£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 (AE5), 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 The iteration process of this example can be considered a model of the density-based method; if the initial value of Δ 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.
In the following part, the same problem is solved with different initial designs to illustrate the initial designdependent 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 Δ ¼ 0: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. 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 Φ 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. 5.2 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.
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. Δ 0 denotes the initial Δ value, dΔ denotes the reduction value of Δ in each iteration, and minΔ denotes the minimum value of Δ. N indicates the total iteration number when the optimization convergences, and MC means the mean compliance of the design.
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 designdependent 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.

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. [59]. 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.

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 Φ 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 Φ 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.