Efficient Local Level Set Method without Reinitialization and Its Appliance to Topology Optimization

The local level set method (LLSM) is higher than the LSMs with global models in computational efficiency, because of the use of narrow-band model. The computational efficiency of the LLSM can be further increased by avoiding the reinitialization procedure by introducing a distance regularized equation (DRE).The numerical stability of the DRE can be ensured by a proposed conditionally stable difference scheme under reverse diffusion constraints. Nevertheless, the proposed method possesses no mechanism to nucleate new holes in the material domain for two-dimensional structures, so that a bidirectional evolutionary algorithm based on discrete level set functions is combined with the LLSM to replace the numerical process of hole nucleation. Numerical examples are given to show high computational efficiency and numerical stability of this algorithm for topology optimization.


Introduction
Topology optimization is a numerical iterative procedure for making an optimal layout of a structure or the best distribution of material in the conceptual design stage [1].The level set method (LSM) is a recently developed approach to topology optimization that uses a flexible implicit description of the material domain [2].The central idea of the LSM is to employ an implicit boundary describing model to parameterize the geometric model, and the boundary of a structure is embedded in a high-dimensional level set function that is called its zero level set [3].The level set-based method is able to not only fundamentally avoid checkerboards and mesh-dependence, but also maintain smooth boundaries and distinct material interfaces during the topological design process [4].Hence, many level set-based methods [5] have been developed for topology optimization since the LSM was first introduced into structure optimization.
With an implicit local level set model, the computational efficiency of the local level set method (LLSM) [6] is much higher than that of the global level set methods, especially for shape optimization.However, the main shortcoming of the conventional LSM is that it possesses no mechanism to nucleate new holes in the material domain for twodimensional structures, resulting in the final design heavily dependent on the initial guess [4].A mechanism named the bubble-method [7] was first proposed to create new holes inside the structures in topology and shape optimization.This idea has been further developed into the mathematical concept of topological derivatives [8].In the shapesensitivity-based level set approaches, topological derivatives are incorporated to indicate the best place for introducing a new hole in a separate step of the optimization process [9] or as an additional term in the Hamilton-Jacobi equation [10].The globally supported radial basis function (RBF) [11] and compactly supported RBF (CSRBF) [12] are typically used to discretize the original time-dependent initial value problem into an interpolation problem.The CSRBF brings about the strictly positive definiteness and sparseness properties of matrices under certain conditions.Hence the CSRBF has generalized the practical applications of RBFs to a larger set of scattered data [12].
In the conventional LSM [3], a reinitialization procedure usually needs to reshape the level set function (LSF) to a signed distance function (SDF) periodically.However, the zero level set may drift away from its initial position by 2 Mathematical Problems in Engineering iteratively solving a classical reinitialization equation [13].To suppress this drift, an interface preserving level set redistancing algorithm is proposed by Sussman and Fatemi [14].Nevertheless, it has been proved that the SDF is not a feasible solution to the H-J equation [15].In practice, it not only raises serious problems as when and how it should be performed, but also affects numerical accuracy in an undesirable way and thus should be avoided as much as possible [16].The need for reinitialization was originally eliminated by introducing a penalty term [17] into a variational level set approach [18].The undesirable boundary effect of the penalty term can be eliminated by taking a distance regularized equation (DRE) instead of this term.Hence, a so-called distance regularized level set evolution (DRLSE) [16] is realized based on the variational approach.As an unnecessary diffusion effect of the DRLSE was found in some locations where the surface is too flat, the DRE was recently modified with a new and balanced formulation to eliminate this effect [19].Although parts of the diffusion rates in the DRE are negative, the numerical stability can still be maintained by incorporating reverse diffusion constraints in the difference schemes of the DRE, as can the reverse diffusion equations with all negative diffusion rates [20].
The aim of this work is to solve the aforementioned numerical issues that still exist in the LLSM for topology optimization of two-dimensional structures.A bidirectional evolutionary algorithm based on the discrete level set functions (DLSFs) is proposed to find a stable topological solution first and then combined with the LLSM to further evolve the local details of the topology and shape of the structure.Transforming the DLSFs into the local level set function of the LLSM is achieved by iteratively solving the DRE.After that, the DRE is incorporated into the LLSM to avoid the reinitialization procedure.A difference scheme under reverse diffusion constraints is formulated for the DRE to improve its numerical stability.Typical examples are given to show the effectiveness of the proposed algorithm in terms of convergence, computational efficiency, and numerical stability.

Local Level Set Method
Using Narrow-Band Model.In the local level set method (LLSM) [6], the local level set equation is defined as where (, ) is defined as the local level set function (LLSF) and   is the normal velocity in normal direction  = ∇/|∇|; the truncation function () is with  0 = { :     ()     < } being a narrow band with the half-band width .The narrow-band model and the corresponding LLSF are described as shown in Figure 1.
It can be seen from Figure 1 that only the LLSF within the narrow-band  0 needs to be updated during each iteration.Hence, the LLSM is higher in computational efficiency than the LSMs based on global level set models.

Bidirectional Evolutionary
where  0 is a predefined constant set as 1 in this study.
The values of elemental densities can be derived from the DLSFs.If the th element belongs to  1 , that is,  ∈  1 , then the element density   = 1; if  ∈  3 , then   =  min , where  min is a small value 0.001; if  ∈  2 , then   ∈ ( min , 1), where   is calculated in terms of the interpolation criterion given in the code manual [21].In the structural model, the rectangular element is split into four triangles first, and the value of the DLSF at the common point of the triangles is then given by the average of the values of the four points.After that, each triangle is examined separately in the same logic.Finally, the elemental density is found to be the average of the contributions of the triangles.
The structural stiffness design has been widely investigated in numerous literatures for topological sensitivity analysis.The standard notion [1] of minimum compliance design problems under a global volume constraint can be mathematically defined as follows: where  is known as the mean compliance, the open set Ω represents all admissible shapes in the design region , () is the nodal displacement vector, and  denotes the global stiffness matrix;   is the volume of an individual element,  * is the prescribed total volume, and  is the number of elements.Similar to the bubble-method [7] and the level set-based optimization methods [4,10], topological derivatives are taken as topological sensitivities in this study.The topological derivative for node  is given by [9]

Zero level set
where   denotes the material elasticity tensor for node ,  is the strain tensor, and the lame constants with the Poisson ratio ] and Young's modulus of solid materials  0 .tr() denotes the trace of a matrix .
Based on an interpolation function proposed by Shepard [22], a filter scheme of mesh independence is proposed to avoid the checkerboard patterns and mesh dependencies.A circular domain Ω  is first defined as the influence region centered round point x with cut-off radius   , and  Ω denotes the number of points located inside the influence domain.The sensitivity filtering using the Shepard method with scattered points is then defined by where   (x) is the Shepard interpolation with the basis denotes the radial distance from point  to   , and if   (x) ≥   , then   (x) is set to zero. is a positive constant and chosen as a onefold mesh size in terms of numerical experiences.
Over the last two decades, many topology description models have been developed for topology optimization of structures, which can roughly be classified into two categories, the material distribution model and the boundary description model [23].Based on the material distribution model, the ESO (Evolutionary Structural Optimization) method has won a great deal of popularity in recent years [24].The bidirectional ESO (BESO) method [25], as an extension of the ESO method, allows efficient material to be added to the structure while the inefficient one is removed simultaneously.So a bidirectional evolutionary algorithm is developed by integrating both the DLSFs and topological derivatives into the optimization criteria of the BESO method [25].Note that the design variables and topological sensitivities in the BESO method are based on the elemental pseudo densities while those in the proposed algorithm are based on the discrete level set functions.
It is assumed that the volume in the th iteration   is known, and  ≥ 0. The target volume  +1 in the next iteration is then updated by with the evolutionary volume ratio ER and the volume limit  * defined in (4).
A parameter AR  is defined as the adding number of nodes in the set   1 divided by the total numbers of nodes and AR  ≤ AR  max , where AR  max is a predefined positive constant.The definition of AR  is different from that of AR in the original BESO method, since the former parameter corresponds to nodal sensitivity while the latter one corresponds to elemental sensitivity.
It is assumed that the DLSF    of node  is known in the th iteration.Then the DLSF  +1  in the next iteration is updated by where the threshold sensitivity numbers  del th and  add th are determined as the number of nodes decreased from the set   1 and that increased from the set   3 , respectively.These thresholds are similar to those based on elemental sensitivities given in the original BESO method.Full details of determining these thresholds are described in [25].
Finally, a stable topological solution is obtained when the following convergence criterion is satisfied: where  is an allowable convergence error with typical values ranging from 0.001 to 0.01.The majority of logical steps of the bidirectional evolutionary algorithm are presented in Figure 2.
If |∇ 0 | ≤ 0.5 is satisfied for all the initial values  0 , the diffusion effect of (10) will make |∇| approach 0. So it loses the ability to regularize |∇| to 1.An improved diffusion rate  2 () with a diffusion function like the following is proposed in [19]: where  is a positive constant and is chosen as fourfold mesh sizes in terms of numerical experiences. 1 = { :     ()     < } is a narrow band with a half-band width .
The diffusion effect can be divided into two parts: the forward diffusion for |∇| ≥ 1 and the backward diffusion for |∇| < 1.It will make |∇| approach one within  1 but zero outside  1 .However, the two parts are balanced within  1 but unbalanced outside  1 so that multiple iterations are required to retain a flatter level set surface outside  1 .In this paper, the diffusion function  2 (|∇|) is further localized by introducing the half-band width  of the narrow-band  0 in LLSM, thereby resulting in an improved diffusion rate  3 () using the diffusion function With the diffusion rate  3 (), the two parts are balanced within  0 without influencing the level set surface outside  0 .

A Conditionally Stable Difference Scheme for DRE.
It is noted that the common difference schemes for the DRE with parts of the negative diffusion rates are incapable of remaining stable during an iterative process, according to the stability definition of the difference equation.In our numerical experiments, || is apt to gradually become divergent along with the process of iterations.To enhance the numerical  stability of the DRE, a difference scheme similar to that of the mean curvature given in [18] It has been verified by our numerical experiments that the evolution of level set can remain bounded stability even after a large number of iterations for solving (15).However, the maximum of || often exceeds the initial value  0 defined by (1), which leads to the level set surface unsmoothed near the edges of the narrow-band  0 , thereby reducing the computational accuracy of (15).Furthermore, multiple iterations are required to find a suitable Courant-Friedrichs-Lewy (CFL) condition to ensure the stability of (15).The issues related to the numerical instability of the DRE can be resolved by imposing reverse diffusion constraints on the difference scheme (see ( 15)), since it has been proved that the constraints can ensure the numerical stability of the diffusion equations with all negative diffusion rates [20].
First, (15) can be subdivided along the direction of  and  into It can be seen that the four diffusion rates in ( 15) satisfy the boundedness |(  )| ≤ .If  ≤ Δ, the reverse diffusion constraints can be defined by If Δ/Δ ≤ 1/4, first substituting (20) into inequalities (21) and then substituting the result into (18), one can obtain a solution using the inequalities min ,=−1,0,1 It can be seen that the absolute values of   +,+ for ,  = −1, 0, 1 are lower than their initial value  0 .That means that all the absolute values || ≤  0 if the CFL conditions are satisfied:

Flow Chart for Difference Schemes to LLSM with DRE.
The procedure for the LLSM with the DRE consists of two main parts, transforming the models of discrete level set functions into the local level set function and solving the difference schemes of the LLSE associated with the DRE.The final DLSFs corresponding to the stable topological solution can be transformed into the LLSF within the initial narrowband  0 by iteratively solving the DRE under reverse diffusion constraints.The LLSE can be solved by difference schemes using the third-order Runge-Kutta (R-K) scheme for temporal discretization and the fifth-order weighted essentially nonoscillatory (WENO) scheme for spatial discretization.The reader is referred to [26] for more numerical details.The logical steps of the two parts can be described by a flow chart given in Figure 3.  Mathematical Problems in Engineering 2.6.Shape Derivatives and Normal Extension Velocities.With the classical level set model [4], the minimum compliance design problem given in (4) can be converted to an unconstrained problem with the Lagrangian method: where the Lagrangian function (, ) is the objective functional. + is the Lagrangian multiplier of the volume constraint.() is the Heaviside function.The mean compliance (, ) is reformulated as For a number of level set-based approaches [4,9,11,12], the steepest descent method is used to ensure the decrease of the objective function by directly setting the normal velocity field   as the negative shape derivative of (, ).For the particular case of a 2-D model of a linear elastic structure, the boundary traction is fixed and remains unchanged, the displacement constraint is fixed, and the body force is set to zero; thus   can be given by The reader is referred to the article [4] for more detailed theoretical proofs.In addition, a bisectioning algorithm is used to find the Lagrangian multiplier  + to guarantee that the volume constraint be exactly satisfied during each iteration.
The normal velocity field can be naturally extended to the whole domain using the so-called "ersatz material" approach, which fills the void areas with a weak material and then the material density is assumed to be piecewise constant in each element and is adequately interpolated in those elements cut by the zero level set (the shape boundary) [4].In the LLSM, the extension velocity field is localized within the narrowband  0 .By iteratively solving the difference scheme for the DRE (see (15)), one can obtain a smooth velocity field in the region near the edges of the narrow-band  0 to improve the computational accuracy of the extension velocity.

Numerical Examples
In this section, two widely researched examples, the cantilever beam and the arch bridge, are presented in the context of structural minimum compliance design to demonstrate the characteristics of the proposed method.Some of the system parameters using the same values are defined as follows.
Young's elasticity modulus for the solid material is  0 = 200 GPa and for weak material is  min = 10 −3 Pa, and the Poisson ratio is 0.3.The volume constraint  * = 0.5 all , where  all is the total volume in the design region .The convergence tolerance  is set as 0.01 in the bidirectional evolutionary algorithm and 0.001 in the LLSM. 4 is the design domain of a cantilever beam with a size 40 mm × 25 mm.The left side of the domain is fixed as the Dirichlet boundary, and a concentrated force  = 100 N is vertically applied at the central point of the right side as a nonhomogeneous Neumann boundary.In the bidirectional evolutionary algorithm, ER = 2%, AR  max = 5%, and   = 2 mm.In the difference schemes for the LLSE, Δ = 0.1, Δ = 0.3, and  = 0.999.In the difference schemes for the DRE,  = 0.5 and  = 0.001.The design domain is discretized with a mesh of 80 × 50 quadrilateral elements.

First Cantilever-Beam Model. Shown in Figure
In the design domain as shown in Figure 4, the initial volume  0 of the solid region Ω is set to  0 =  all .The structural topologies corresponding to the zero level set and related level set surfaces are shown in Figures 5 and 6, respectively.The convergence histories of the mean compliance and volume fraction are depicted in Figure 7.The result in Figure 5(e) stands for a stable topological solution obtained from the bidirectional evolutionary algorithm.Topological results given by this algorithm are characterized by a smooth boundary attributed to the structural model described by the DLSFs.By comparing Figures 6(e) and 6(f), the sharp level set surface corresponding to the DLSFs has been successfully converted into a smooth one related to a local level set function by iteratively solving the DRE at the initial stage of the LLSM.Then the shape of the boundary is further improved by iteratively solving the LLSE.In addition, all the absolute values of the level set function are less than the initial value  0 .Therefore it verifies the effectiveness of reverse diffusion constraints on the numerical stability of the difference scheme for the DRE.

Second Cantilever-Beam
Model.This study has also investigated the influence of different initial models of structure on the final design.Figure 8 depicts the design domain with a size 4.0 mm × 2.5 mm.The left side of the domain is fixed and a concentrated force  = 1 N is vertically applied at the bottom of the right side.All the parameters but   = 0.2 mm and AR  max = 1% remain unchanged as those of the first cantilever-beam model.Figure 9 shows two cases of the initial configurations with full materials and the least but essential materials and their topologies during the process of optimization.The two final designs are made with the same topology and almost similar shape of the structure, which shows the complexity of the final topology is not changing appreciably with different initial structures.Therefore, the numerical process of the bidirectional evolutionary algorithm can be used to replace the numerical process of hole nucleation in the LLSM to avoid the final design heavily dependent on the initial guess.
Figure 10 shows the topological topologies for almost   = 0.2 mm using several mesh sizes.It can be seen that the optimal topology does not depend on the discretization in terms of layout and number of bars.side, and the sum of the pressure is 5 N.In accordance with the same definitions of the cantilever-beam parameters, the parameters are set to ER = 3%, AR  max = 5%,   = 0.1 mm, Δ = 0.02, Δ = 0.3,  = 0.99,  = 0.001, and  = 0.0167.The design domain is discretized with a mesh of 120 × 72 quadrilateral elements.
This example focuses on the new characteristic of the proposed algorithm for improving the convergence of the bidirectional evolutionary algorithm using the LLSM.The structural topologies corresponding to the zero level set are depicted in Figure 12.Note that it starts from the initial model with the volume  = 0.5 all and remains unchanged, so as to maintain the stability of the evolution process of the object function.The evolutionary histories of the objective and the volume constraint starting from the initial models are plotted in Figure 13.A design of the structure shown in Figure 12(b) corresponding to a stable topological solution is also the final design of the topology (not shape) made merely using the bidirectional evolutionary algorithm in our numerical experiments.The subsequent topologies given in Figures 12(c)-12(f) show that the LLSM can further optimize the topology of the structure to improve convergence.Moreover, the LLSM can also improve the shape of the boundary to obtain a smoother shape design till it reaches the convergence tolerance in the 38th iteration.It is worth noticing that the final topology obtained by the LLSM is just a local optimal design because of the use of the steepest descent method.Despite the optimal solution of this arch-bridge model obtained by an element-wise ESO method [27], in this case the optimized topology obtained

Initial
Step 20 Step 40 Step 43 (a) The case starts from the full-material initial configurations

Initial
Step 20 Step 60 Step 76 (b) The case starts from the least-material initial configurations   LLSM, a bidirectional evolutionary algorithm is combined with the LLSM.This proposed algorithm has been used successfully in topology optimization of two-dimensional (2-D) structures, and it is easy to be extended to 3-D structures.The main features of this algorithm unknown to the conventional LSMs and the LLSM can be summarized as follows: (a) The discrete level set functions can be efficiently transformed into the local level set function by iteratively solving the distance regularized equation (DRE).
(b) The DRE can be used instead of the reinitialization equation to further increase the computational efficiency of the LLSM.
(c) A conditionally stable difference scheme under reverse diffusion constraints is formulated to ensure the numerical stability of the DRE.
(d) If the stable topological solutions of the bidirectional evolutionary algorithm are inconsistent, the LLSM can achieve at least the consistent local optimal solution for the different cases of initial models.High computational efficiency and numerical stability of the proposed algorithm have been verified by three typical numerical examples.

Figure 1 :
Figure 1: Narrow-band model and local level set function in LLSM.

j
Carry out the finite element analysis; calculate the topological derivative of each node the initial DLSFs  0 Define the elemental set S and for the next iteration V k+1 Update the DLSFs  k+1 j according to the Stop iteration and obtain number k = k + 1

Figure 2 :
Figure 2: Flow chart depicting logical steps of the bidirectional evolutionary algorithm.

Figure 3 :
Figure 3: Flow chart depicting logical steps of the LLSM.

3. 3 .Figure 4 :Figure 5 :
Figure 4: Design domain of the first cantilever beam and its boundary conditions.

Figure 9 :Figure 10 :Figure 11 :
Figure 9: Topologies of zero level set in the two cases of the second cantilever beam.

Figure 13 :
Figure 13: Convergent histories of the objective function and the constraint.

Figure 14 :
Figure 14: Topologies of zero level set for the initial model, and the final models obtained by the bidirectional evolutionary algorithm and the LLSM, respectively: (a) case 1,  = 0.4 all ; (b) case 2,  = 0.6 all ; (c) case 3,  = 0.8 all .
Algorithm with Discrete Level Set Functions.A two-dimensional structural model is built in the work region  ⊂  2 .And the set  =  1 ∪  2 ∪  3 represents the finite elements in the domain .It can be divided into three parts,  1 that consists of the solid elements with a full-material density;  2 that covers the elements with intermediate material densities;  3 that involves the void elements with a weak-material density.Accordingly, the nodal sets corresponding to the elemental sets  1 ,  2 , and  3 is developed and described as ) +1/2,   +  , −  (  ) −1/2,   −  , Δ +  (  ) ,+1/2   +  , −  (  ) ,−1/2  2 .