Abstract
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.
1. 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 two-dimensional 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 shape-sensitivity-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 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 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.
2. Optimization Algorithm
2.1. Local Level Set Method Using Narrow-Band Model
In the local level set method (LLSM) [6], the local level set equation is defined aswhere is defined as the local level set function (LLSF) and is the normal velocity in normal direction ; the truncation function is with 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 needs to be updated during each iteration. Hence, the LLSM is higher in computational efficiency than the LSMs based on global level set models.
2.2. Bidirectional Evolutionary Algorithm with Discrete Level Set Functions
A two-dimensional structural model is built in the work region . And the set represents the finite elements in the domain . It can be divided into three parts, that consists of the solid elements with a full-material density; that covers the elements with intermediate material densities; that involves the void elements with a weak-material density. Accordingly, the nodal sets corresponding to the elemental sets , and are defined as , , and , respectively. If it is assumed that , consist of all the nodes within the region , then a discrete level set function (DLSF) for node can be defined as where is a predefined constant set as 1 in this study.
The values of elemental densities can be derived from the DLSFs. If the element belongs to , that is, , then the element density ; if , then , where is a small value 0.001; if , then , 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]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 . 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 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 bywhere is the Shepard interpolation with the basis function , in which denotes the radial distance from point to , and if , then 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 iteration is known, and . The target volume in the next iteration is then updated bywith the evolutionary volume ratio and the volume limit defined in (4).
A parameter is defined as the adding number of nodes in the set divided by the total numbers of nodes and , where is a predefined positive constant. The definition of is different from that of 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 iteration. Then the DLSF in the next iteration is updated by where the threshold sensitivity numbers and are determined as the number of nodes decreased from the set and that increased from the set , 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.
2.3. Distance Regularized Equation (DRE) and Its Improvement
In the distance regularized level set evolution (DRLSE) [16], the DRE can retain the signed distance feature at least within the narrow-band region near boundaries without reinitialization, whose formula is expressed in the standard form of the diffusion equation aswith the diffusion rate , where the diffusion function is set to with .
In the original DRLSE, the energy density was defined as which is a double-well potential function because there are two minimum points of at and . So the diffusion function is given by
It is easy to verify the boundedness of the diffusion rate and . It can be seen from (12) that, for , is positive, and will decrease and approach 1; for , is negative, and will increase and approach 1; for , is positive, and will decrease and approach 0.
If is satisfied for all the initial values , the diffusion effect of (10) will make approach 0. So it loses the ability to regularize to 1. An improved diffusion rate 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. is a narrow band with a half-band width .
The diffusion effect can be divided into two parts: the forward diffusion for and the backward diffusion for . It will make approach one within but zero outside . However, the two parts are balanced within but unbalanced outside so that multiple iterations are required to retain a flatter level set surface outside . In this paper, the diffusion function is further localized by introducing the half-band width of the narrow-band in LLSM, thereby resulting in an improved diffusion rate using the diffusion function
With the diffusion rate , the two parts are balanced within without influencing the level set surface outside .
2.4. 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] is developed and described aswhereand , in which the difference schemes of are given by
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 defined by (1), which leads to the level set surface unsmoothed near the edges of the narrow-band , 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 intowith and .
Then the four flow functions are defined aswhere denotes the change from to in one time step and in the direction, and the definitions of , , and are similar to that of . The lowest and highest limit values of these flow functions are defined as
It can be seen that the four diffusion rates in (15) satisfy the boundedness . If , the reverse diffusion constraints can be defined by
If , first substituting (20) into inequalities (21) and then substituting the result into (18), one can obtain a solution using the inequalities
It can be seen that the absolute values of for are lower than their initial value . That means that all the absolute values if the CFL conditions are satisfied:
2.5. 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 narrow-band 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.
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 narrow-band . 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 to improve the computational accuracy of the extension velocity.
3. 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 and for weak material is , and the Poisson ratio is 0.3. The volume constraint , where 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.
3.1. First Cantilever-Beam Model
Shown in Figure 4 is the design domain of a cantilever beam with a size . The left side of the domain is fixed as the Dirichlet boundary, and a concentrated force is vertically applied at the central point of the right side as a nonhomogeneous Neumann boundary. In the bidirectional evolutionary algorithm, , , and . In the difference schemes for the LLSE, , , and . In the difference schemes for the DRE, and . The design domain is discretized with a mesh of quadrilateral elements.
In the design domain as shown in Figure 4, the initial volume of the solid region is set to . 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 . Therefore it verifies the effectiveness of reverse diffusion constraints on the numerical stability of the difference scheme for the DRE.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
3.2. 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 . The left side of the domain is fixed and a concentrated force is vertically applied at the bottom of the right side. All the parameters but and 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.
(a) The case starts from the full-material initial configurations
(b) The case starts from the least-material initial configurations
Figure 10 shows the topological topologies for almost 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.
(a)
(b)
(c)
(d)
3.3. Arch-Bridge Model
The design domain of an arch-bridge model with a size is shown in Figure 11. Both the bottom corners of the domain are the fixed support. A uniform static pressure is vertically applied on the upper 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 , , , , , , , and . The design domain is discretized with a mesh of 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 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 by the proposed bidirectional evolutionary algorithm is not reasonable compared with the final optimal solution shown in Figure 12(f).
(a)
(b)
(c)
(d)
(e)
(f)
Starting from different initial models, the final models obtained by the bidirectional evolutionary algorithm and the LLSM, respectively, are shown in Figure 14.
(a)
(b)
(c)
It can be seen from Figures 12 and 14 that the final topologies obtained by the bidirectional evolutionary algorithm are inconsistent. In contrast, the final optimized results found using the LLSM subsequently are of the same topology and similar shape. Although the local optimal solution of the arch-bridge model can also be obtained by using the ESO/BESO methods with elemental variables, numerical instabilities and zigzag boundaries can result in these elemental variables-based methods. Hence, nodal variables are needed to take the place of the elemental variables in these methods. The combined algorithm with the bidirectional evolutionary algorithm and the LLSM can also resolve this problem and achieve at least the consistent local optimal solution for the different cases of initial models.
4. Conclusions
The LLSM is intended to remarkably increase the computational efficiency of the conventional LSMs using global models. To overcome the issue of hole nucleation of the 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of the paper.
Acknowledgment
The financial support from National Natural Science Foundation of China (no. 51278218) is gratefully acknowledged.