Efficient Local Level Set Algorithm Combined with Bidirectional Evolutionary Criterion Using Discrete Variables and its Appliance to Topology and Shape Optimization

Many level set methods (LSMs) have been used for topology and shape optimization. An efficient topology optimization algorithm is proposed by combining a new bi-directional evolutionary criterion (BEC) with the local LSM (LLSM). The BEC is established using multiple discrete variables (DVs) and introduced into the bi-directional evolutionary structural optimization (BESO) method so that the DV-BESO method obtained can be applied to topology optimization. The LLSM is used instead of the global LSM due to its high efficiency. In the LLSM, the computational efficiency can be further improved by solving a distance regularized equation without the re-initialization procedure. Two benchmarks show the high efficiency of this proposed combined algorithm for topology and shape optimization.


Introduction
It is widely known that many level set-based methods [1] have been developed for topology and shape optimization since Sethian et al. [2] first introduced the LSM into structure optimization. By using a local level set (LLS) model, the LLS method (LLSM) [3] is with much higher convergence efficiency than the LSMs with global models. For the conventional LSMs and the LLSM, new holes cannot be generated in the material region far away from structural boundaries. For the global LSMs, topological derivatives [4,5] can be introduced to improve the hole-nucleation capability, but cannot be applied to generate holes inside the LLS model. Another challenge is to increase the computational efficiency of the LSMs by avoiding the procedure of reinitializing the signed distance function (SDF) [6], which needs multiple iterations for preserving the shape of the SDF. To overcome this challenge, a distance regularized LSM (DRLSM) [7] can maintain the signed distance feature at least in the local region across the interface, and whose boundary effect is further eliminated by improving the distance regularized equation (DRE) [8].
Besides the LSMs, the evolutionary structural optimization (ESO) method is also popular for topology optimization, whose most obvious characteristic is only two discrete values of design variables such as 0 and 1. Then the bi-directional ESO (BESO) method [9] is developed to improve the convergence of the ESO method. Both the methods are classified as the hard-kill methods [10] because they remove completely the void elements by setting the density of these elements as 0. However, these hard-kill methods have received considerable criticism for their drawbacks such as lack of rigorous proof of optimality criteria, non-optimal solution, and so on [11]. Two types of BESO methods can be verified to overcome these drawbacks. The first is the so-called soft-kill BESO (SBESO) method [12]. The SBESO model is similar to the SIMP (Solid Isotropic Material with Penalization) model. There are still only two discrete values of elemental densities while the smaller one is a minimum constant (eg. 10 -3 ), not 0. The second is the BESO method with multiple discrete variables (DVs) [13], which is called the DV-BESO method for short. The DV-BESO model is similar to the SBESO model except for the different discrete values of elemental densities and the corresponding bi-directional evolutionary criteria (BEC). In the DV-BESO method, the BEC is not yet sufficiently researched for topology optimization [14].
The aim of this study is to propose an efficient algorithm for topology and shape optimization, by combining the DV-BESO method with the LLSM. In the DV-BESO method, multiple DVs are taken as design variables and an improved BEC is established to update the DVs efficiently. For increasing efficiency, we adopt the LLSM and take the DRE to avoid the reinitialization procedure. Two types of benchmark examples are implemented by the proposed combined algorithm for topology and shape optimization.

Optimization problem and its analytical sensitivities
The structural stiffness design has been widely investigated in numerous literature for topological sensitivity analysis. The standard notion [15] of minimum compliance design problems under a global volume constraint can be mathematically defined as follows: where C is known as the mean compliance, u is the nodal displacement vector, and K denotes the global stiffness matrix; V is the structural volume fraction, V* is the prescribed total volume fraction, ρe is the discrete elemental density with the minmum contant ρmin.
Ghabraie [16] For the special cases that ρe = ρ1 and ρe = ρrmax, only one of equation ( Then the filter scheme proposed in the SBESO method [12] can be extended to ∆ + and ∆ − for the numerical stability, in which a filter radius rmin is defined as the cutoff radius of a circular domain centered at element e. The readers can refer to the equation (16) and (17) of the literature [12] for details.

Bi-directional evolutionary criterion of DV-BESO method
Using multiple DVs, the BEC developed in the DV-BESO method is different from the BEC used in the SBESO method. In the SBESO method, the void element with ρe = ρmin can be converted into the solid element with ρe = 1 in one iteration. In contrast, the BEC of the DV-BESO method is more complicated to update multiple DVs. we take five DVs for example, set ρe = [ρ1, ρ2, ρ3, ρ4, ρ5] = [ρmin, 0.25, 0.5, 0.75, 1], so get ρ1 < ρ2 < ρ3 < ρ4 < ρ5.
In the kth iteration, the current structural volume fraction is set as Vk, two of the evolutionary volume ratios er and arr are defined for the cases of volume decreasing and volume increasing, respectively. In the kth iteration, the structural volume is given by Then we find the upper limit of the increasing number Nadmax = floor(NallVkARmax) and that of the decreasing number Ndemax = floor(NallVkDelmax), respectively, where the function floor(*) denotes rounding down the value in brackets, the sum of elements is Nall, as well as the upper limit of the increasing ratio is ARmax and that of the decreasing ratio are Delmax.
The implementation procedure of the BES is described as follows.
Firstly, is the number of the added elements with the discrete densities ρe ={ρ1, ρ2, ρ3}, and = { , }. We select elements with the larger sensitivity number ∆ + , and increase the corresponding densities by one level, ie. ρe is increased from ρ1 to ρ2, and so on.
Secondly, all the elements are sorted by their sensitivity numbers ∆ − . The threshold sensitivity number ∆ ℎ − is found by the bisection method, and we decrease the densities of the elements whose sensitivity numbers satisfy ∆ − < ∆ ℎ − by one level, ie. ρe is reduced from ρ5 to ρ4, and so on. The volume frction V' is obtained by the bisection method until it satisfies |V' -Vk+1| ≤ ε, where ε is a predefined minimal constant (eg. 10 -3 ).
(2) Vk ≤ V* Firstly, is the number of the reduced elements with ρe ={ρ3, ρ4, ρ5}, and = { , }. We select elements with the smaller sensitivity number ∆ − , and decrease the corresponding densities by one level.
Secondly, all the elements are sorted by ∆ + . ∆ ℎ + is found by the bisection method, and we increase the densities of the elements whose sensitivity numbers satisfy ∆ + > ∆ ℎ + by one level, and V  is determined by the bisection method until it satisfies |V' -Vk+1| ≤ ε.
The BEC usually drives discrete densities from the two endpoint values ρ1 and ρ5 towards the middle values ρ2 ~ ρ4 in each iteration, so that there are too many elements with ρe = ρ2 ~ ρ4 in the final topology. Therefore, a heuristic operator of enhanced contrast scale [17] has been developed to suppress grey-scale boundaries : The parameter c1 is increased gradully, the enhanced effect of the operator T1 also becomes larger. If c1 →∞, T1 will degenerate into a threshold function and there are only two cases for the DVs, ie. ρe = ρmin and Finally, a convergence criterion is required to stop the iterative process. By referring to the SBESO method, we define the convergence error τ, set the convergence criterion Different from the SBESO method, two values of τ need to define as τ0 = 0.001 and τ1 = 0.00001. We firstly set τ = τ0, and introduce the operator of enhanced contrast scale T1(ρe) to drive the elemental densities ρe gradually towards 0/1 distribution when the convergence criterion is satisfied, then terminate the iterative process until equation (7) with τ = τ1 is satisfied .

Procedure of DV-BESO algorithm
The evolutionary iteration procedure of the DV-BESO method for minimizing the mean compliance while satisfying the volume limit is organized as follows: Step 1: Discrete the whole design domain using a finite element mesh.
Step 3: Determine the target volume for the next design Vk+1 according to equation (5).
Step 5: Update multiple discrete variables using the proposed bi-directional evolutionary criterion.
Step 6: If the convergent criterion defined in equation (7) with τ = τ0 is satisfied, update ρe using the operator of enhanced contrast scale T1. Otherwise, go to Step 3.
Step 7: If the objective volume is achieved and the convergent criterion defined in equation (7) with τ = τ1 is satisfied, terminate the iterative process. Otherwise, go to Step 3. Finally, a stable topological solution can be obtained by the DV-BESO method.

Local level set method with distance regularized equation
After the optimization process of the DV-BESO method, the LLSM can be applied to further improve the local topologies and structural boundaries. However, there are great differences between the design variables of the above methods. It needs a transformation process from the DVs crooesponding to the stable topological solution into the initial local level set (LLS) function (LLSF) crooesponding to the LLS model.

Initialization of discrete level set function
Firstly, the nodal density ρn of node n is obtained by an interpolation function proposed by Shepard [18].  (8) where We(xn) is the Shepard interpolation with the basis function Di(xn), in which di(xn) = ||xn − xn|| denotes the distance from the center of element i to node n and it satisfies Di(xn) = 0 for di(xn) ≥ rw where rw is the cutoff radius of a circular domain Ωw centered at node n, Nn represents the number of nodes located inside Ωw, the positive constant c is chosen as 1.5 times of the mesh size. The initial values of discrete LSF (DLSF) 0 is set as −c0 and c0, and then we define 0 at node n as where ϕ(x, t) is defined as the LLSF, Vn is the normal velocity in normal direction n = ϕ/|ϕ|; the truncation function C(ϕ) is with T0 = {x:||ϕ(x)|<γ}is a narrow band with the half-band width γ. The narrow band model and the corresponding LLSF are described as follows. It can be seen from figure 1 that only the LLSF within the narrow band T0 needs to be updated during each iteration. Hence, the LLSM is higher in computational efficiency than the LSMs based on global level set models.
The classical level set method [19] has derived the shape derivatives of the minimum compliance design problem. The readers can refer to the literature [19] for specific details.
In the DRLSM [11], the DRE can retain the signed distance feature at least within the narrow-band region near boundaries without reinitialization. In our earlier study [20], we have discussed the boundary effect of the DRE and its improved strategy [12], construct the final form of the diffusion equation described as with the diffusion rate α(ϕ) = μdp(|ϕ|), where the diffusion function dp(|ϕ|) is ( ) In addition, we define the initial version of the diffusion rate α0(ϕ) = μdp0(|ϕ|) with ( ) ( ) ( ) For the LLSM, lots of difference schemes to equation (10) have been well-developed. For the LLSM, lots of difference schemes to equation (10) have been well-developed. According to the report [21], we adopt the third-order Runge-Kutta (R-K) scheme for temporal discretization and the fifth-order weighted essentially non-oscillatory (WENO) scheme for spatial discretization.
For the DRE, conventional difference schemes to equation (12) are theoretically unstable because some of the diffusion rates are negative. Our earlier study [20] has studied the improved difference schemes to the DRE. The readers can refer to the equation (15) to (23) of the literature [20] for details.

Procedure of trasformation process and local level set method
The transformation process from the DVs e  into the LLSF  is organized as follows: Step 1: Transforming the DVs ρe into nodal densities ρn by introducing equation (8).
The LLSF obtained is limited within the range of T0, and the improved difference schemes to the DRE can ensure its numerical stability.
The evolutionary iteration procedure of the LLSM is organized as follows: Step 1: Calculate the current mean compliance and volume fraction V.
Step 2: Carry out FEA, calculate the shape derivatives of the mean compliance.
Step 3: Solve iteratively the LLSE (equation (10)) with the DRE (equation (12)) until the pre-defined conditions are satisfied. During the step, the total volume fraction can remain unchanged, ie. V = V*.

MBB beam using 80 50
 mesh Firstly, we establish a famous benchmark model-MBB beam, and the right half of 40mm×25mm domain (with 1 mm thickness) is discretized by a 80×50 four-node quadrilateral element mesh; a sliding support along the x-direction is applied to the left boundary, a simple support is applied in the bottom right corner, an external force P = 100N is applied to the top left corner; the material has Young's modulus E = 200GPa. Secondly, several essential parameters need be defined : ρmin = 0.001, the filter radius for the sensitivity number rmin = 3, the filter radius for the Shepard interpolation function rw = 1.5; the evolutionary ratio er = arr = ARmax = Delmax, er0 is the initial value of er, er0 = 0.1∆ρe, er = er0/bb by starting from the 200th iteration and update bb = bb + bbadd in every 20 iterations, and bb is initialized as 2.5 and bbadd = 2.5; for the operator of enhanced contrast scale, c1 is initialized as 1 by starting from the 200th iteration and update c1 = c1 + 0.05 in every 20 iterations; ∆ = 0.3 and γ = 0.999 are set in the LLSE; the target volume fraction V* = 0.3. Figure 2 describes the topological results of the proposed algorithm combined with the DV-BESO method and the LLSM in three stages. The objective is to minimize the mean compliance while the structural volume is limited. For the discrete elemental densities ρe, the number of discrete values Ns = 11, the variation of discrete density ∆ρe = 1 according to the equation ∆ρe = 1/(Ns -1). Firstly, the DV-BESO method was used to find a topological stable solution described by the discrete density model given in figure  2(a). Furthermore, the discrete density model was transformed into the initial model of the LLSF ϕ. The corresponding interpolation model of elemental densities is shown in figure 2(b). The existence of the narrow band T0 can be verified from the shaded part of figure 2(b). Finally, the LLSM was combined to find the topology and shape optimization solution described the 2-dimension(2D) LLS model given in figure 2(c). The 3D surface model of the LLSF is described in figure 2(d). It can be verified from figure 2(d) that the LLSF ϕ is only active within T0 while ϕ = 0 outside T0; the DRE can keep the smooth distribution of ϕ at least near the structural boundaries without reinitialization; the absolute value of ϕ does not larger than c0 defined in equation (9), thus verifying the numerical stability of difference schemes of the DRE improved by our earlier work [20]. The convergent histories for the mean compliance and the volume fraction are depicted in figure 3. It can be seen from the curve of mean compliance that there are two sections with the oscillation characteristic. The first section represents the convergent process of the DV-BESO method, and the convergent history is gradually becoming more stable until it turns into a straight line. The second section represents the process of transform (PoT) from ρe to ϕ ( ρe →ϕ ). By introducing the operator of enhanced contrast scale, the elemental densities corresponding to the stable topological solution tends to be the 0/1 distribution. Firstly, ρe is transformed into the nodal densities ρn by solving the Shepard interpolation function, resulting that the mean compliance is greatly reduced because of the generation of lots of elements with intermediate densities. Secondly, ρn is transformed into ϕ. The PoT ρn →ϕ also can generate the elements as above while these elements is much less than those corresponding to the PoT ρe →ρn, thus the mean compliance corresponding to ϕ gradually approaches to that corresponding to ρe.
As shown in figure 3, the minimum compliance corresponding to the LLSM is less than that corresponding to the DV-BESO method. It is verified that the LLSM can find the better topology and shape optimization solution than that obtain by the DV-BESO method.

MBB beam using 160 80
 mesh In the case of multiple values of ∆ρe, the second example studies comparatively the above two methods on the final topologies and mean compliance as well as the influence of the PoT ρe →ϕ on the optimization process and the topological results. To show more details of structural topology for the optimization results, a modified MBB beam with 40mm×20mm domain needs be discretized by a 160×80 four-node quadrilateral element mesh, other parameters are the same as the last example. The topological results  As shown in figure 4, almost all the minimum compliance can be found on the PoT ρe →ϕ in each case of ∆ρe. However, the final topology corresponding to the PoT ρe →ϕ cannot be considered as the optimal solution. There are two main reasons given as follows: (1) the evolutionary direction of the PoT ρe →ϕ is not necessarily along the optimization direction due to its heuristic characteristic; (2) the PoT ρe →ϕ can generate the elements with intermediate densities, which may result in the smaller mean compliance than these based on the DVs ρe or the LLSF ϕ. Therefore, the final topologies corresponding to the LLSM can be still taken as the topological optimal solutions.
It is necessary to illustrate several beneficial effects of the PoT ρe →ϕ. Firstly, the PoT ρe →ϕ can enhance the possibility of jumping out of the original final topology obtained by the DV-BESO method. For a small value of V*, there will inevitably be a few thin components, which are difficult to manufacture. These thin components can be quickly deleted by introducing the transformation ρe →ϕ, but not all of them will be eliminated, otherwise, it will be affected that the convergent stability of the optimization process, and the total volume will not maintain balance. Secondly, the PoT ρe →ϕ helps the LLSM to find the better topologies. In most cases, the LLSM can only improve structural topologies indirectly by the evolution of structural boundaries, and its ability to topology modification is not as good as that of the DV-BESO method; the PoT ρe →ϕ is used to transform the elemental density-based model with the 0/1 distribution to the LLS model within the narrow band T0, which makes it easier to improve local topologies and converge to a new topology by the LLSM. Thirdly, the PoT ρe →ϕ can reduce greatly the mean compliance without or only a little modification for structural topologies, so the final compliance corresponding to the LLSM is generally less than that corresponding to the DV-BESO method.

Conclusions
The LLSM using the local model is intended to remarkably increase the computational efficiency of the conventional LSMs using global models. However, the LLSM cannot generate new holes inside the structural domain. In contrast, it has been verified that the effectiveness of the BESO method for topology optimization. A combined algorithm is proposed with the BESO method and the LLSM. The main features of this algorithm unknown to the conventional BESO method and the LLSM can be summarized as follows: a) a bi-directional evolutionary criterion with discrete variables (DVs) is proposed instead of the optimization criterion of the BESO method, the DV-BESO method obtained can be used to find a stable topological solution; b) a heuristic strategy is proposed to transform the DVs of the DV-BESO method into the local level set function (LLSF) of the LLSM by using the proposed interpolation function and the distance regularized equation (DRE); c) the LLSM is combined with the DV-BESO method to further improve local topologies and the shape of structural boundaries, and its capability to topolog8.y optimization can be enhanced by the transformation from the DVs into the LLSF; d) the DRE can be used instead of the reinitialization procedure to further increase the computational efficiency of the LLSM.
The high computational efficiency and numerical stability of the proposed algorithm have been verified by two benchmark examples.