Parametric structural shape & topology optimization with a variational distance-regularized level set method

Accepted Manuscript Parametric structural shape & topology optimization with a variational distance-regularized level set method Long Jiang, Shikui Chen PII: DOI: Reference:

S0045-7825(16)31406-2 http://dx.doi.org/10.1016/j.cma.2017.03.044 CMA 11401

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date: 20 October 2016 Revised date: 29 March 2017 Accepted date: 31 March 2017 Please cite this article as: L. Jiang, S. Chen, Parametric structural shape & topology optimization with a variational distance-regularized level set method, Comput. Methods Appl. Mech. Engrg. (2017), http://dx.doi.org/10.1016/j.cma.2017.03.044 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Parametric Structural Shape & Topology Optimization with a Variational Distance-Regularized Level Set Method Long Jiang and Shikui Chen 1 Computational Modeling, Analysis and Design Optimization Research Laboratory Department of Mechanical Engineering State University of New York Stony Brook, N.Y, 11794 ABSTRACT The signed distance function (SDF) gives the shortest distance from a given point to the boundary, and the sign indicates whether this point is inside or outside the closed boundary or enclosed region. The SDF property is highly preferred in classical level set methods to maintain the numerical stability during the topology optimization process and provide a metric for the distance-based interpolation of different material properties. In conventional level set methods, a common way of achieving a level set function with the signed distance property is to periodically implement the so-called reinitialization scheme by solving an additional Hamilton-Jacobi partial differential equation. However, such reinitialization scheme is implemented outside the optimization loop with the optimization process suspended, which may shift the optimization result and bring convergence issues. In this paper, a double-well potential functional is employed for distance regularization inside the topology optimization loop, which can enforce the signed distance property of the level set function in a narrow band along the design boundaries while keeping the level set function flat in the rest area of the computational domain. The radial basis function (RBF) based parameterization technique is combined with mathematical programming to improve the performance of the proposed distanceregularized topology optimization method in handling problems with non-convex objective functions and multiple constraints. The flatness of the level set function in the material region also enables easy creation of new holes to the design in the topology optimization process. Both 2D and 3D benchmark examples are employed to demonstrate the validity of the proposed method. Keywords: Parametric level set method; Distance regularization; Topology optimization

1

Address all correspondence to this author. The State University of New York at Stony Brook. Department of Mechanical Engineering. Tel: +1 (631) 632-2309 Fax: +1 (631) 632-8544. Email address: [email protected] (S. Chen)

1

1. INTRODUCTION Topology optimization is a numerical scheme that can find the optimal material distribution in a design domain, subject to specified boundary conditions by minimizing (or maximizing) the objective function and satisfying one or multiple design constraints [1]. Since Bendsøe and Kikuchi proposed the homogenization-based topology optimization approach in 1988 [2], the research on structural shape and topology optimization has experienced a boost in the last decades [3]. As a variant of the homogenization method, Solid Isotropic Material with Penalization (SIMP) method [4] is widely employed in many applications due to its simplicity and computational efficiency. In SIMP, an artificial density, varying between 0 and 1, is mapped to the material property by a power law, and the optimal (or optimized) structure layout can be achieved by finding the distribution of the artificial density field, with a pre-specified objective function minimized or maximized. Despite its easy-to-implement advantage, the SIMP method often cannot result in a design with sharp boundaries but a design with intermediate densities [5]. The intermediate density, also termed the gray area or grayscale elements, may blur the boundaries and smear the details of the optimization result, making it difficult for design verification and fabrication. The intermediate density result in the conventional density-based method is unavoidable because of the original mathematical implementation of the method itself [6]. Many researchers proposed different schemes to handle this issue and generate clear final designs. Those schemes include, but not limited to regularized intermediate density control [7], perimeter control [8, 9], post-processing schemes [10-12], design feature control [13, 14], projection schemes [15, 16] and generating final designs via pre-defined clear-boundary elements [17], and so on. However, the downside of introducing extra schemes is that they could make the result deviate away from the optimal design, and the performance of the original design may degenerate after post-processing. Unlike the homogenization or the SIMP method, the level set approach is a “boundarybased” approach which directly evolves the design boundaries during the topology optimization and guarantees a final design with clear boundaries [18]. Level set methods were originally proposed by Osher and Sethian [19] for tracking the evolution of moving fluid fronts, where the merging and splitting of fluids lead to topological changes. Level set methods implicitly represent the design boundaries as the zero level set of the onehigher-dimension level set function [18]. During the topology optimization process, a design velocity field, usually calculated by shape sensitivity analysis [20], is introduced to the Hamilton-Jacobi partial differential equation (PDE) to evolve the level set function. With level set methods, the boundary locations can be captured at sub-element accuracy, and other geometric information, such as normal vectors or curvatures, are also embedded in the level set model. Osher and Santosa [21] introduced the shape gradient of the objective functional into the level set model, setting up a link between the shape gradient and the design velocity field. This work is further advanced by Allaire [20, 22] and Wang [23, 24]. Allaire et al. derived the shape sensitivity of compliance and the geometric advantage by employing the adjoint variable method. Following the material derivative method in continuum mechanics. Wang et al. [23] identified a link between the design velocity filed in the level set methods and the general structural sensitivity analysis. To have more control over the final design with better manufacturability, Mei et 2

al. [25] introduced the feature-based topology optimization for the structural design approach. A real value function was employed to represent the design structure boundaries in an implicit way. A recent contribution from Guo et al. [26], extended the conventional level set method by imposing the constraints on the extreme value of the sign-distanced level set function, so the feature of the final design can be controlled. Compared with the earlier work of Guest et al. [27] about the design feature controlling, Guo et al.’s approach eliminated the need for introducing the multiple nonlinear feature control constraints, thereby gaining higher computational efficiency. A similar approach was made by Allaire et al. [28] by introducing geometric manufacturing constraints in the structural optimization. The min/max local structure thickness and the minimum members’ distance constraint were considered, and they were all calculated via the Heaviside function of the design shape. By using a specific integral criterion, the satisfaction of the proposed constraints was monitored. With the powerful tool of the level set methods, other types of constraints can be considered when doing shape and topology optimization in different areas. Some recent investigations are done on the casting constraint proposed by Wang et al. [29], the geometric constraint in architectural design proposed by Dapogny et al. [30], the arbitrary geometric constraint in the bone structure design proposed by Deng et al. [31] and so on. Despite its advantages over SIMP or homogenization method, deficiencies of classical level set approach cannot be ignored. To be particular, the following three basic issues need to be considered when the classical level-set approach is applied to topology optimization [32]. First, the level set function may become too flat or too steep during the evolution, since for a given zero level set, the corresponding level set function is not necessarily to be unique. The flatness or the steepness will lead to undesired numerical stability issues. In conventional level set methods, this issue is alleviated by periodically reinitializing the level set function to be a signed distance function, which essentially regulates the function value at a given point to be the shortest distance to the boundary. The sign of the distance marks whether the point is inside or outside the closed design boundary. This scheme is called the “reinitialization” [33, 34] and it usually requires solving another Hamilton-Jacobi PDE [35] or directly computing the distance based on the geometric boundary information [34]. However, this conventional reinitialization scheme is implemented outside the topology optimization process, which needs to suspend the optimization and therefore is less efficient. Even so, when and where to impose the reinitialization is not certainly known yet [36]. The periodical reinitialization of the level set function will smear out the newly formed design details, i.e. holes inside the material region, and will decrease the topological flexibility of the level set scheme. To help create new holes during the level set evolution, Dunning et al. [37] utilized a secondary level set function to facilitate new hole insertion. However, a better hole creation scheme should be developed to let the level set function itself generate the final design result in a natural way. Second, to meet the Courant-Friedrichs-Lewy (CFL) necessary condition for stability and convergence [38], the limit of the time step in the explicit time integration scheme should be small enough, so in each iteration the boundary will evolve no more than one grid interval. In another word, to achieve an accurate result, the mesh size should be 3

small enough. This, in turn, will require more iterations to reach the final design and much more computing time in each iteration. Last but not the least, to evolve the level set function, the design velocity field must be extended to the whole design domain or at least to a narrow band along the boundary. However, the design velocity field obtained from shape sensitivity analysis is nominally defined only at the design boundary. In classical level set approach, this usually brings the need of solving an additional Hamilton-Jacobi PDE [39] to extend the velocity field from the boundary to the remaining area of the computational domain, which is not desired in the optimization process. To overcome the limitations of the conventional level set approach, many researchers have proposed different strategies to perfect the level set approach. To increase the computational efficiency and design flexibility, Guo et al. [40] proposed a structural topology optimization approach based on Moving Morphable Components (MMC). In this approach, the design elements were multiple components explicitly defined by several control parameters, such as the position, thickness, or angle, etc. The control parameters were updated during the optimization to change the design accordingly. Later, to adjust the geometric complexity of the final design, Yamada et al. [41] introduced an additional fictitious energy term to the original objective function of structural topology optimization. Similarly, introducing a distance-suppression energy term to force the level set function to be a signed distance function was explored by Zhu et al. [42]. Otomori et al. [43] extended Yamada’s work to the design of compliant mechanisms by introducing mathematical programming. Helped by the mathematical programming, the solving of the Hamilton-Jacobi Partial Differential Equation (PDE) can be avoided and introducing multiple constraints become easier. Moreover, mathematical programming can better improve the computational efficiency compared with the conventional level set methods. Many researchers have made remarkable contributions to increase the robustness and efficiency of level set approach in solving non-convex optimization problems with multiple constraints. Wang and others [44, 45] proposed to parameterize the level set function with radial basis functions (RBFs), which transformed the original HamiltonJacobi PDE to a system of ordinary differential equations (ODEs) to reduce the computational burden. Shapiro et al. [46] proposed a similar method by representing the level set function using B-splines and use R-functions to support desired parametric changes. Later, Qian et al. [47] parameterized the structural design with Non-Uniform Rational B-Splines (NURBS) for better numerical accuracy and computational efficiency in geometric conversion from the design model to the analysis model. Both positions and weights of the NURBS were design variables in shape optimization to increase the flexibility and robustness of shape representation. A more thorough investigation of topology optimization via B-Splines is summarized in Qian’s later work [48] with more examples. Recently, Zhang et al. [49] not only incorporated the B-splines to represent the smooth and clear design boundary but also introduced the Moving Morphable Components (MMC) to increase the adaptively of the design result. Being an explicit topology optimization methods, MMC can provide direct control of the geometric features, such as the radius of an individual hole, the length, thickness, or rotation angle of a structural component, which is an obvious advantage over implicit approaches. The RBF or B-splines parameterization scheme provided some promising features for 4

topology optimization process [50], such as better robustness in optimization, no need of reinitialization, and more flexibility in generating new holes. Once the level set function is parameterized with RBFs or B-splines, gradient-based mathematical programming optimizers, such as CONLIN [51], Method of Moving Asymptotes (MMA), or GCMMA [52, 53], can carry out the optimization. It is worth noting the parametric topology optimization method no longer regularizes the level set function to be a signed distance function, although a signed distance level set function is always preferred when implementing finite element analysis on a fixed mesh. The reason a signed distance type level set function is preferred is discussed as follows. Even though the level set representation can always separate the “material region” from the “void region”, it is not always the case with the material property interpolation. In the level set based material interpolation scheme, whether a specific area of the design domain is a material domain or not is usually determined by the Heaviside function value of the level set function. Mathematically, the width of the transition zone for the Heaviside function cannot be smaller than an element size, so there is always an area along the “material region” that is a “gray area”. If the level set function over the transition zone is flat, then the gray area after material interpolation may be relatively large. However, when the level set function is relatively steep, the gray area can be minimized. But it may take a long time for those steep surfaces atop zero level set to move downwards and generate new holes. That is why the level set function should be maintained in the desired shape during the optimization process. Otomori et al. [43] noticed this issue and alleviated the gray area by introducing a relatively small transition width for the Heaviside function of the level set function. In this paper, the authors propose to use a new variational level set method which can intrinsically maintain the shape of the level set function throughout the optimization process. A double-well potential functional originally proposed by Li et al. [54] for image segmentation is employed here to achieve distance-regularized level set evolution. By minimizing the double-well potential functional, the norm of the gradient of the level set function can be regularized to be 1 (a signed distance shape) within the narrow band along the design boundary while kept being 0 (a flat surface) elsewhere. Compared with the single-well energy functional utilized in the work of Zhu et al. [42], where the norm of the gradient of the level set function is only kept being 1, the double-well energy function gives an extra feature to the level set function of being flat in the non-transition zone, making it easier for hole generation during the optimization. Moreover, the signed distance property along the boundary will ensure an accurate material interpolation. The desired distance-regularized level set evolution can increase the stability, robustness and topological flexibility of the optimization. The nonlinear potential term for distance regularization can either be added to the objective functional as a penalty or be treated as a constraint of the structural optimization problem setting. When the potential term is added as part of the objective functional, a new weighted objective functional can be formed up, so the level set function can maintain the distance-regularized shape during the optimization process. Treating the potential term as a constraint can also enforce the desired distanceregularized property of the level set function in each step of the optimization process. In Zhu’s work [42], the single-well energy functional is only introduced in the objective 5

function, and the update of the level set function is done by solving a generalized Hamilton-Jacobi PDE, which can be less computational efficient compared with mathematical programming. In this paper, the structural topology optimization problem is solved via an RBF-based parametric level set method (PLSM). The sensitivity of the mean compliance, the distance-regularized potential, and the constraints are derived via shape sensitivity analysis. With the RBF-based parameterization of the level set function, the free form boundary is driven by the mathematical programming optimizer while the topological flexibility in conventional PDE-based level set methods is retained. The double-well potential energy functional is firstly minimized along with the original objective functional with a penalty and later constrained to be zero in the optimization process to demonstrate its robustness. By combining the distance-regularization energy functional with the parametric level set method and mathematical programming, the advantage of clear boundary representation of the level set methods can be kept, while the unwanted features of conventional level set methods such as the reinitialization scheme, the constraint of CFL condition and the velocity extension can be eliminated. The remaining paper is arranged as follows: In section 2, the concept of the conventional level set methods is introduced. The distance regularization energy with the conventional parametric level set evolution is detailed in section 3. Section 4 performs the structural topology optimization problem with the conventional parametric scheme. In section 5, the regularization term with the parameterization scheme is combined with the structural topology optimization problem to show the distance regularization effect during the optimization process. This proposed technique is employed as a penalty in the original objective functional in section 5.1 and as a constraint as well in section 5.2 to verify the performance. The conclusions are drawn in section 6. 2. CONVENTIONAL LEVEL SET TOPOLOGY OPTIMIZATION

METHODS

FOR

STRUCTURAL

The general idea of level set methods for the structural topology optimization problem is to represent the desired structural boundaries as the zero contour (in 2D) or isosurface (in 3D) of a level set function with one higher dimension [18]. In this paper, we use the symbol D to represent the design domain, Ω to stand for the material region, and Γ to denote its boundaries. The level set representation can be defined as Φ ( x, t ) > 0, ( x ∈ Ω )  ( x, t ) 0, ( x ∈ Γ ) , Φ= Φ ( x, t ) < 0 , ( x ∈ D \ Ω ) 

(1)

where x is an arbitrary point in the design domain D ; t is the pseudo time [20, 24] of the dynamic shape optimization process.

6

Figure 1. The level set function and its corresponding zero level set

0 on both sides with respect to By differentiating the boundary expression of Φ ( x, t ) = a pseudo time t , a Hamilton-Jacobi PDE can be achieved as equation (2), ∂Φ − Vn ∇Φ = 0 , ∂t

(2)

where Vn is the normal velocity field calculated through sensitivity analysis. The level set based structural optimization was firstly proposed by Sethian and Wiegmann [55] using the level set boundary representation mentioned above. Based on it, the minimum mean compliance structural topology optimization can be formulated as follows: Minimize := J (u , Φ ) Φ

∫ j ( u ) H ( Φ )d Ω D

s. t.

∫

Volume =

D

u |∂= u0 Du

H ( Φ )= d Ω V0 ,

(3)

∀v ∈ U ,

a ( u , v, Φ= ) l ( v, Φ ) In equation (3), the objective functional J ( u , Φ ) represents the mean compliance of the whole structure. The strain energy density j ( u , Ω ) can be defined as Eijkl ε ij ( u ) ε kl ( u ) where Eijkl stands for the elasticity tensor and ε ij ( u ) is the strain tensor. u is denoting the displacement field in the kinematically admissible displacement fields U and u0 is the prescribed displacement. The final volume of the optimization is targeted to be V0 . In

equation (3), a ( u , v, Φ ) is the bilinear form of the energy and l ( v, Φ ) the linear form of the load, which are formulated as follows: = a ( u , v, Φ ) = l ( v, Φ )

∫

D

∫

D

Eijkl ε ij ( u )ε kl ( v ) H ( Φ ) d Ω

pvH ( Φ ) d Ω + ∫ τ vδ ( Φ ) ∇Φ d Ω

,

(4)

D

where p is the body force and τ is the boundary traction. 7

In equation (3) and (4), H ( Φ ) represents the Heaviside step function, which is set to be zero in the void area and one in the material region. In practice, the Heaviside function is smoothed within a narrow band in the transition area, since a step function is not differentiable at the point where H ( Φ ) is discontinuous. A common selection of the band width is between two or three times of the spatial mesh size. One of the smoothed Heaviside functions [24] can be formulated as:

δ, Φ < −∆    Φ Φ3  1 + δ  H= , −∆ ≤ Φ ≤ ∆ . 0.75 (1 − δ )  − 3  + δ (Φ)  ∆ 3∆  2     1, Φ>∆ 

(5)

Here δ is a small number while ∆ is half of the transition width. The original Heaviside function and its smoothed version are illustrated below:

Figure 2. Smoothed Heaviside function (dashed line) and original Heaviside function (solid line)

By applying the Lagrange multiplier method for this minimum compliance problem, the Vn can be expressed as

Vn =−λ + Eijkl ε ij ( u ) ε kl ( v )

(6)

For the details of the derivation of the sensitivity, please be referred to [20, 24] 3. DISTANCE-REGULARIZED EVOLUTION PARAMETRIC LEVEL SET METHOD

LEVEL-SET

WITH

A

3.1 The distance regularization energy functional The conventional level set needs the reinitialization scheme to regularize the level set function to be a signed distance function during the optimization process. In this paper, the reinitialization scheme is replaced by employing a distance-regularizing energy functional from [54] to regularize the level set evolution during the optimization process. A general the distance-regularizing energy functional R can be formulated as follows, 8

R  ∫ P ( ∇Φ )d Ω , Ω

(7)

where P ( ∇Φ ) stands for the regularization energy potential; Φ denotes the level set function and ∇Φ represents the norm of ∇Φ . To maintain the signed distance property of the level set function, ∇Φ = 1 needs to be satisfied. Based on this need, if the regularization term, which is a function of ∇Φ , is to be minimized along with the original objective functional, the global minimum of the regularization term should be reached when ∇Φ = 1 . In [54], Li et al. proposed two types of regularization potential with the desired characteristics, which are named as single-well potential and double-well potential. They are denoted as P1 and P2 in equation (8) and (9), respectively: 2 1 ∇Φ − 1) , ( 2  1  1 − cos ( 2π ∇Φ )  , ∇Φ < 1 2    ( 2π ) P2 =  . 2 1  ∇Φ > 1 ( ∇Φ − 1) ,  2 The corresponding plots of the two potentials are provided in Figure 3.

= P1

(8)

(9)

Figure 3. (a) The single-well potential P1 and its derivative P1 ' . (b) The double-well potential

P2 and its derivative P2 ' .

The derivatives of the potentials are P1 ' =

∂P1 ∂P2 and P2 ' = , respectively. ∂ ∇Φ ∂ ∇Φ

As illustrated in Figure 3 (b), the double-well potential P2 has two global minimum points at ∇Φ =1 and ∇Φ = 0 . With the potential minimized, the level set function can

9

be forced either to be a signed distance function ( ∇Φ =1 ) within the narrow band along the design boundary or to be flat surfaces ( ∇Φ =0 ) elsewhere. This is the so-called distance-regularized level set function. With the level set function in the narrow band regularized to be a signed distance function, the mapping from the level set function to the structural material properties will be more accurate. At the same time, the flatness of the level set function will make easier nucleation of new holes during the optimization process. The generation of new holes will add extra flexibility and robustness to the level set evolution. This is a property that the conventional level set does not possess since the reinitialization scheme for numerical stability control will smear out the newly generated holes when the level set function is forced to be a signed-distance function. This shows the necessity of introducing the distance regularization term to the optimization process. It can be seen from Figure 3 (b) that within the range from 0 to 1, the double-well energy functional has the local maximum at ∇Φ = 0.5 and the functional is symmetrical within this range. The local maximum works as the selecting criterion to determine whether the value of the function ∇Φ will be forced to be 0 or to 1. A further modification of the double-well potential energy function can be made to meet different selecting requirements. 3.2 RBF-based parameterization In this section, the radial basis function (RBF) is introduced to parameterize the level set function. The RBF is formulated based on its center at the point with a radiallysymmetric property. Generally, an RBF can be defined as the following form [56]:

(

)

= φi φ x − xic .

(10)

Here  denotes the Euclidean norm and xic is the location of the selected center. There are a number of different RBFs, including Gaussian, Multi Quadric (MQ), Inverse Multi Quadric (IMQ), Compactly Supported RBFs (CSRBF) and so on. The readers can refer to [44, 57, 58] for details about different RBFs. In this paper, the CSRBF with C2 continuity is chosen due to its strictly positive definiteness, sparseness, and smoothness [56]. Given the normalized distance r , the CSRBF with C2 continuity can be formulated as below:

φ (r ) = max{0, (1 − r ) } ( 4r + 1) 4

where r =

(x − x ) +( y − y ) c 2 i

c 2 i

.

(11)

d min

Here the normalized distance r is also called the support radius, is defined by dividing the distance between the node ( x, y ) and the center node ( xic , yic ) , by a searching range of d min . This CSRBF is illustrated in Figure 4 with a support radius of 0.8 in a 2-by-2 square domain. 10

Figure 4. The CSRBF with C2 continuity

For the selection of RBF under different scenarios and the identification of the parameters for different RBFs, the readers are referred to [59] for details. With a specified RBF: φ ( x ) , the level set function Φ can be decomposed into a sum of multiplication between RBFs on the nodes and their corresponding weights: Φ ( x, t ) = ∑ α i ( t ) φi ( x ) . n

i =1

(12)

Here the weight of the ith node, α i , is also known as the expansion coefficient and φi is the selected RBF on this node. The CSRBF is a spatial function, and the expansion coefficient α i is only time dependent. Substitute equation (12) into equation (2), the original Hamilton-Jacobi PDE is parameterized into the ODE form: ∑ α i ( t ) φi ( x ) − Vn ∇Φ = 0 . n

i =1

(13)

When performing the CSRBF based parameterization scheme, a proper support radius for the CSRBF must be selected. It should not be too small to span enough neighboring nodes, nor too big in order not to increase the computational cost [60]. But theoretically, a relatively big support radius leads to design with a smoother boundary. Based on the level set representation, an unconstrained topology optimization design problem is set up as follows considering a design-independent load [44]: = Min L ( u , Φ ) Φ( x)

1 Eijkl ε ij ( u ) ε kl ( u ) H ( Φ ) d Ω + λ ∫ H ( Φ ) d Ω , D 2 ∫D

where λ is a positive Lagrange multiplier and

∫

D

(14)

H ( Φ ) d Ω represents the volume of the

design structure. Moreover, the topology optimization design problem is simplified to be a more computationally efficient unconstrained problem setting, without an explicit final volume target. If a volume target is imposed, the value of the Lagrange multiplier must be updated throughout the optimization process, as detailed in [61]. Following [20, 24], we can get the shape derivative of the above equation as dL = dΩ

∫

Γ

λ − Eijkl ε ij ( u ) ε kl ( u ) Vn ds .

(15) 11

Following the steepest descent method, we can construct the normal velocity Vn as equation (16),

= Vn Eijkl ε ij ( u ) ε kl ( u ) − λ .

(16)

The above normal velocity field is utilized in the conventional level set method [24]. With the parameterization scheme, the Vn in equation (13) can be formulated as: Vn =

1 n ∑ α i ( t ) φi ( x ) . ∇Φ i =1

(17)

By substituting equation (17) into equation (15) and rearranging the position of the time derived expansion coefficient, the equation can be re-written as: dL n 1 = ∑ α i ( t ) ∫ λ − Eijkl ε ij ( u ) ε kl ( u )  φi ( x ) ds . Γ dt i =1 ∇Φ

(18)

By applying the chain rule, equation (18) can be also written as: n ∂V dL n ∂J = ∑ αi + λ ∑ αi . = i 1 ∂α ∂α i dt i 1 = i

(19)

By comparing the corresponding parts of equation (18) and (19), the sensitivity can be formulated as follows [62]: ∂J 1 = − ∫ Eijkl ε ij ( u ) ε kl ( u ) 1,..., n φi ( x ) ds , i = Γ ∂α i ∇Φ . (20) ∂V 1 = ∫= i 1,..., n φ ( x ) ds , Γ ∇Φ i ∂α i The above boundary integration can be transformed into a domain integration utilizing equation (21),

ds= δ ( Φ ) ∇Φ d Ω ,

(21)

where δ () is the Dirac delta function. 4. STRUCTURAL TOPOLOGY OPTIMIZATION WITH CONVENTIONAL PARAMETRIC LEVEL SET METHOD In equation (13), the conventional Hamilton-Jacobi PDE is parameterized into an ordinary differential equation (ODE). Changing the values of the expansion coefficient α i will change the shape of the level set function Φ . In the parametric level set method, the expansion coefficients are selected to be the design variables, to be updated through the topology optimization process. Due to the smoothness that the selected CSRBF with C2 continuity possesses, the mature gradient based optimization solver, Method of Moving Asymptotes (MMA) [52], is utilized in the optimization process to generate accurate and convincing results. 12

A minimum compliance of a cantilever beam problem is presented here to demonstrate the conventional parametric level set method. A 2-by-1 rectangular domain is fixed on the left edge with a point load F = 1 applied at the center of the right edge. The mesh size is set to be 0.02 with an 8 times mesh size support radius of the CSRBF. The volume target is aimed at 0.5. To separate the material domain from the void domain, the values of the nominal Young’s modulus for the elastic and dummy materials are set to be 100,000 and 1 respectively. The Poisson’s ratios of both elastic and dummy materials are set to be 0.3. The figure below shows the problem boundary conditions.

Figure 5. Problem settings for the cantilever beam benchmark example

With the sensitivity analysis results in section 3, the optimization result generated by the parametric level set method (PLSM) is shown below:

Figure 6. The convergence history for PLSM

13

Figure 7. Topology optimization process of a cantilever beam using conventional parametric level set method (PLSM). Top row: level set functions; and bottom row: initial (left), intermediate (middle) and final (right) designs

As seen from the result, the objective function is minimized with the volume constraint well satisfied. However, it can be seen that after some iterations the level set function is subject to a drastic fluctuation throughout the computational domain. This numerical unsteadiness should be eased in order not to accumulate any potential numerical errors. Combining the parametric level set methods with the above distance regularization energy can ensure a relatively steady level set evolution through the optimization process. The details are introduced in the next section. 5. STRUCTURAL TOPOLOGY OPTIMIZATION WITH REGULARIZED PARAMETRIC LEVEL SET METHOD

DISTANCE-

5.1 Distance regularization energy utilized as an objective function penalty In this section the distance regularization energy functional is combined with the parametric level set method for structural topology optimization, aiming to achieve the distance-regularized level set evolution throughout the optimization process. The level set function is parameterized using the CSRBF with C2 continuity, and the expansion coefficients are updated as the design variables during the optimization process. The distance-regularized effect is illustrated below in a geometric manner. As shown in Figure 8 (a), the characters “NY” is built in a BMP black-white image. This image is read into the program next to be treated as the initial level set function, as shown in Figure 8 (b). By minimizing the distance regularization energy, the transition narrow band of the level set function is regularized into the signed distance shape. The level set function at places far from the transition narrow band is preserved as flat surfaces, as shown in Figure 8 (c). The shape of the level set function illustrated in Figure 8 (c) is the preferred distance regularized level set function expected to be achieved with the 14

proposed scheme. With this distance-regularized level set function, the numerical advantages mentioned in the introduction can be secured. The readers can compare Figure 8 with Figure 1 to have a better understanding of the distance-regularized level set function.

Figure 8. The visualization of the distance regularization effect. (a) The initial BMP image. (b) The initial level set function (c) Distance-regularized level set function

With the parametric level set method, the structural optimization problem can be solved with MMA or other mathematical programming optimizers. The optimization formulation is similar to equation (3) except the regularization term R with its penalty w in the objective function: Minimize= : J ( u,α ) α

∫ j ( u ) H ( Φ )dΩ + wR D

sub. to :

∫

D

(22)

H ( Φ )d Ω =Vconst ,

a ( u,v= u0 ) L ( v ) u |∂= Du

∀v ∈ U

The sensitivity analysis is given as follows. By applying the chain rule, the derivative of the distance regularization term regarding α i can be formulated as

∂R ∂R ∂Φ = ⋅ = − div  d p ( ∇Φ ) ∇Φ  ⋅ φi ( x ) . ∂α i ∂Φ ∂α i

(23)

Here the φi (x) is the selected RBF (CSRBF with C2 continuity) and div is the divergence operator. d p is a function defined as [54]:

15

d p (s) 

p′( s ) s

(24)

The p′( s ) is the first derivative of energy functional p ( s ) . In this paper the p ( s ) is the double-well potential P2 shown in equation (9). Then equation (23) can be plugged into the sensitivity in equation (20) with a penalty w to be used to update the parametric level set function: ∂J 1 ∂R = − ∫ Eijkl ε ij ( u ) ε kl ( u ) φi ( x ) ds + w , i= 1,..., n Γ ∂α i ∇Φ ∂α i . (25) ∂V 1 = ∫= φ ( x ) ds , i 1,..., n Γ ∇Φ i ∂α i The rectangular design domain in section 4 is used again. The vertical force applied at the center of the right edge is amplified to be 10. The left edge of the design domain is fixed. The design domain is discretized by the mesh size of 0.02 with the support radius of 8 times of the mesh size. The values of the nominal Young’s Modulus are set to be 100,000 and 1 for the elastic and dummy material. With the Poisson’s ratio to be 0.3 for the elastic and dummy materials and the final volume targeted at 0.45, the following result can be got.

Figure 9. The convergence history for DRPLSM

16

Figure 10. Topology optimization process of a cantilever beam using distance-regularized parametric level set method (DRPLSM). Top row: level set functions; and bottom row: initial (left), intermediate (middle) and final (right) designs

The weighted objective function is minimized with the satisfaction of the volume ratio target at 0.45. No huge fluctuation can be seen from the level set evolution anymore, and the level set function is maintained to be the distance-regularized shape at the end of the optimization. From the third row of Figure 9, the regularization energy is maintained at a relatively low level, proving that the distance-regularized level set function is achieved during the optimization process. Since the level set function is updated via the change of the expansion coefficient without manually imposed reinitialization, theoretically new design details, i.e. new holes can be created freely during the optimization process. The creation of new holes can increase the robustness of the design process, making it preferable for topology optimization. To demonstrate the capability of the hole generation, the previous problem formulation is kept while only to start from an initial design with one hole. New holes are expected to be created during the optimization process. The result is shown below.

17

Figure 11. The convergence history for DRPLSM starting from one hole

Figure 12. Demonstration of the flexibility in hole generation of the proposed DRPLSM. Top row: level set functions; and bottom row: initial (left), intermediate (middle) and final (right) designs

For the intermediate design, it can be noticed that new holes are generated in the optimization process. With the minimization of the combined objective function and the satisfaction of the volume constraint, the final design matches the design result that starts with multiple holes, as in the previous example. 18

Following the similar settings, the distance-regularized parametric level set topology optimization is employed to a 2D MBB beam benchmark example. Due to the symmetry, the design domain is halved to a 2-by-1 rectangle, half of the original 4-by-1 design domain, to increase the computational efficiency. The material properties and the mesh size are kept the same as the previous distance regularized cantilever beam example. The force applied at the upper center segment is set to be F = 10 downwards. The volume target is set to be 0.45. The details of the boundary conditions are shown in Figure 13.

Figure 13. The design domain and boundary conditions of the MBB beam

The optimization result is shown as follows.

Figure 14. The convergence history for MBB beam of DRPLSM

19

Figure 15. Topology optimization process of a 2D MBB beam with the proposed DRPLSM. Top row: level set functions; and bottom row: initial (left), intermediate (middle) and final (right) designs

A Michell-type structure benchmark example with fixed-fixed supports is demonstrated next with multiple loads. A 2-by-1 design domain is discretized with the mesh size of 0.02. Three vertical loads, as shown in Figure 16, with F1 = 10 and F2 = F3 = 5 are applied at the lower edge. The material properties are kept the same as the previous distance-regularized cantilever beam example while the volume target is set to be 0.3. The boundary conditions and the optimization results are shown below:

Figure 16. The design domain and boundary conditions of the Michell-type structure

20

Figure 17. The convergence history for Michell-type structure of DRPLSM

Figure 18. Topology optimization process of a 2D Michell-type structure with the proposed DRPLSM. Top row: level set functions; and bottom row: initial (left), intermediate (middle) and final (right) designs

The proposed method is also extended to a 3D structural topology optimization problem to demonstrate its robustness. The 3D example is a cubic box domain with the edge length equals to 1. A 1000 force is applied at the middle of the lower right edge with the left corners and two edge centers fixed. The domain is discretized with the mesh size of 0.025 with the volume target at 0.1. Due to the symmetry, only half the box (separated 21

along the red plane with slash lines) is utilized, and the result is mirrored back to illustrate the final design.

Figure 19. The design domain and boundary conditions of the 3D cantilever beam

Figure20. The optimized 3D cantilever beam and its regularization details using DRPLSM

In Figure 20, the green isosurface represents the zero level set. The dark iso-surface inside the red circle represents an inner offset of the zero level set iso-surface. It can be seen that the zero level set iso-surface and its offset are parallel with one the other, proving the regularization effect in 3D cases. Besides different boundary conditions, some other parameter settings in the regularized parametric level set methods might influence the result. One of them would be the support radius of the CSRBF. A too small support radius cannot span effectively through the neighbors of the CSRBFs, and a too large radius would increase the computational cost, but with a better smoothing effect. In the previous listed 3D example, the support 22

radius is selected to be 8 times of the mesh size to achieve the balance between numerical accuracy and computational cost. A 25% raise in the support radius of the CSRBF for the 3D cantilever beam example, namely 10 times of the mesh size, is optimized again for comparison purpose.

Figure 21. The 3D cantilever beam with different support radius of DRPLSM

As seen on the right side of the above figure, the new optimization results in a similar final design with the same regularization effect mentioned previously. However, the design generated by a bigger support radius of the CSRBF tends to have smoother boundaries, with fewer branch details. This is because the result is smoothed by a larger overlapping area between the neighbors of adjacent CSRBFs 5.2 Distance regularization energy works as a constraint From the results in Section 5.1, it can be seen that putting the external energy into the objective function is a feasible solution to achieve a distance-regularized level set function due to the specified global minimum of the double-well potential functional. However, selecting the weight ratio w requires a parametric study to have a balance between the regularization effect and the original objective functional. According to equation (22) and (25), the penalty will not only affect the objective function value but also affect the sensitivity of the whole optimization process. This means if the penalty becomes too big, the level set function will be “over-regularized”, preventing generation of the result. If the penalty is given too small, then the distance-regularized effect will not be obvious. Moreover, as the optimization process goes on, the level set function is supposed to change every iteration. Putting the distance regularization energy functional term in the objective function cannot guarantee that the energy term is at its minimum in every single iteration. This means with the distance regularization energy working as an objective functional penalty, a distance-regularized level set function evolution may not be guaranteed in every iteration. To maintain a distance-regularized level set evolution through the whole optimization process, without selecting the penalty value, the regularization term can be introduced as a constraint. Thanks to the advantages mention in the introduction for the RBF 23

parameterization and the gradient based solver MMA, the nonlinear constraint can be easily added to the optimization process. The problem can be formulated as:

Minimize : = J ( u, α ) α

∫ j ( u ) H ( Φ )dΩ D

sub. to :

∫

D

H ( Φ )d Ω =V0 ,

(26)

R = 0, a ( u , v= u0 ) L ( v ) u |∂= Du

∀v ∈ U

One thing must be pointed out is that MMA is not very suitable for handling equality constraints. However, in order to let the final volume ratio hit the target, the authors introduced two opposite inequality constraints together inside MMA. For example, if the volume constraint of the design is targeted at V0 , as ∫ H ( Φ )d Ω =Vconst , then in the D

optimization optimizer MMA, the volume constraint is handled as

∫

D

H ( Φ )d Ω ≤ Vconst

− ∫ H ( Φ )d Ω ≤ −Vconst

.

(27)

D

In this case, the volume of the design is “forced” to reach the volume target. The distance regularization constraint is also handled like this. By treating the regularization term as a constraint, the original 2D cantilever beam example from section 5.1 is utilized again with the volume ratio targeted at 0.45. The sensitivity for the extra constraint can be easily taken from equation (23). The results are shown as follows.

24

Figure 22. The convergence history of the 2D cantilever beam with regularization energy in constraint

Figure 23. The topology optimization process of a cantilever beam using distance-regularized parametric level set method (DRPLSM) with distance regularization energy in the constraint. Top row: level set functions; and bottom row: initial (left), intermediate (middle) and final (right) designs

It can be seen from the convergence plot in Figure 22 that the double-well energy is maintained relatively close to its original value throughout the whole process, ensuring the distance-regularized feature of the level set function and eliminating the need for selecting a proper weight ratio in the objective function. 25

6. CONCLUSION In this paper, a variational distance-regularized parametric level set method is proposed by combining the conventional parametric level set method with a distance regularization energy functional. The proposed method is tested with 2D and 3D structural topology optimization examples. Results show that the desired distanceregularized level set evolution is achieved and the unwanted features in conventional level set methods, like reinitialization, CFL condition, and others are eliminated. Besides that, new holes can be created during the optimization process, compared with the conventional level set methods in which the hole generation can be difficult. The combination of the parametric level set method with the distance regularization term is implemented both as part of the objective functional and as a constraint to show the robustness of the proposed method. This newly proposed method can eliminate the numerical fluctuation existed in the conventional non-regularized parametric level set while avoiding the grayscale issue found in the previous research. The distanceregularized parametric level set method is expected to address more complex structural topology optimization problem, which is also our goal for the next step.

Acknowledgments The authors want to thank the three anonymous reviewers for their helpful and constructive comments that greatly contributed to improving the quality of the paper. This research was supported by the National Science Foundation (CMMI1462270), the unrestricted grant from Ford Research & Advanced Engineering, and the start-up funds from the State University of New York at Stony Brook.

26

REFERENCE [1] M.P. Bendsøe, O. Sigmund, Topology Optimization-Theory, Methods, and Applications, (2003). [2] M.P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Computer Methods in Applied Mechanics and Engineering, 71 (1988) 197-224. [3] O. Sigmund, K. Maute, Topology optimization approaches, Structural and Multidisciplinary Optimization, 48 (2013) 1031-1055. [4] M.P. Bendsøe, Optimal shape design as a material distribution problem, Structural Optimization, 1 (1989) 193-202. [5] T. Bruns, A reevaluation of the SIMP method with filtering and an alternative formulation for solid–void topology optimization, Structural and Multidisciplinary Optimization, 30 (2005) 428-436. [6] M.P. Bendsøe, O. Sigmund, Material interpolation schemes in topology optimization, Archive of Applied Mechanics, 69 (1999) 635-654. [7] T. Borrvall, J. Petersson, Topology optimization using regularized intermediate density control, Computer Methods in Applied Mechanics and Engineering, 190 (2001) 4911-4928. [8] J. Petersson, Some convergence results in perimeter-controlled topology optimization, Computer Methods in Applied Mechanics and Engineering, 171 (1999) 123-140. [9] P. Fernandes, J.M. Guedes, H. Rodrigues, Topology optimization of threedimensional linear elastic structures with a constraint on “perimeter”, Computers & Structures, 73 (1999) 583-594. [10] O. Sigmund, J. Petersson, Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima, Structural Optimization, 16 (1998) 68-75. [11] M.-H. Hsu, Y.-L. Hsu, Interpreting three-dimensional structural topology optimization results, Computers & Structures, 83 (2005) 327-337. [12] A. Koguchi, N. Kikuchi, A surface reconstruction algorithm for topology optimization, Engineering with Computers, 22 (2006) 1-10. [13] J.K. Guest, J.H. Prévost, T. Belytschko, Achieving minimum length scale in topology optimization using nodal design variables and projection functions, International Journal for Numerical Methods in Engineering, 61 (2004) 238-254. [14] T.A. Poulsen, A new scheme for imposing a minimum length scale in topology optimization, International Journal for Numerical Methods in Engineering, 57 (2003) 741-760. [15] F. Wang, B.S. Lazarov, O. Sigmund, On projection methods, convergence and robust formulations in topology optimization, Structural and Multidisciplinary Optimization, 43 (2011) 767-784. [16] J.K. Guest, A. Asadpoure, S.-H. Ha, Eliminating beta-continuation from heaviside projection and density filter algorithms, Structural and Multidisciplinary Optimization, 44 (2011) 443-453. [17] J. Norato, B. Bell, D. Tortorelli, A geometry projection method for continuum-based topology optimization with discrete elements, Computer Methods in Applied Mechanics and Engineering, 293 (2015) 306-327. 27

[18] J.A. Sethian, Theory, algorithms, and applications of level set methods for propagating interfaces, Acta Numerica, 5 (1996) 309-395. [19] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, Journal of Computational Physics, 79 (1988) 12-49. [20] G. Allaire, F. Jouve, A.-M. Toader, Structural optimization using sensitivity analysis and a level-set method, Journal of Computational Physics, 194 (2004) 363-393. [21] S.J. Osher, F. Santosa, Level set methods for optimization problems involving geometry and constraints: I. Frequencies of a two-density inhomogeneous drum, Journal of Computational Physics, 171 (2001) 272-288. [22] G. Allaire, F. Jouve, A.-M. Toader, A level-set method for shape optimization, Comptes Rendus Mathematique, 334 (2002) 1125-1130. [23] M.Y. Wang, X. Wang, PDE-driven level sets, shape sensitivity and curvature flow for structural topology optimization, Computer Modeling in Engineering and Sciences, 6 (2004) 373-396. [24] M.Y. Wang, X. Wang, D. Guo, A level set method for structural topology optimization, Computer methods in applied mechanics and engineering, 192 (2003) 227246. [25] Y. Mei, X. Wang, G. Cheng, A feature-based topological optimization for structure design, Advances in Engineering Software, 39 (2008) 71-87. [26] X. Guo, W. Zhang, W. Zhong, Explicit feature control in structural topology optimization via level set method, Computer Methods in Applied Mechanics and Engineering, 272 (2014) 354-378. [27] J.K. Guest, Imposing maximum length scale in topology optimization, Structural and Multidisciplinary Optimization, 37 (2009) 463-473. [28] G. Allaire, F. Jouve, G. Michailidis, Thickness control in structural optimization via a level set method, Structural and Multidisciplinary Optimization, 53 (2016) 1349-1382. [29] Y. Wang, Z. Kang, Structural shape and topology optimization of cast parts using level set method, International Journal for Numerical Methods in Engineering, (2017). [30] C. Dapogny, A. Faure, G. Michailidis, G. Allaire, A. Couvelas, R. Estevez, Geometric constraints for shape and topology optimization in architectural design, (2016). [31] X. Deng, Y. Wang, J. Yan, T. Liu, S. Wang, Topology optimization of total femur structure: Application of parameterized level set method under geometric constraints, Journal of Mechanical Design, 138 (2016) 011402. [32] Z. Luo, L. Tong, M.Y. Wang, S. Wang, Shape and topology optimization of compliant mechanisms using a parameterization level set method, Journal of Computational Physics, 227 (2007) 680-705. [33] S. Osher, R. Fedkiw, Level set methods and dynamic implicit surfaces, Springer Science & Business Media, 2006. [34] S. Yamasaki, S. Nishiwaki, T. Yamada, K. Izui, M. Yoshimura, A structural optimization method based on the level set method using a new geometry‐based re‐ initialization scheme, International journal for numerical methods in engineering, 83 (2010) 1580-1624. [35] D. Peng, B. Merriman, S. Osher, H. Zhao, M. Kang, A PDE-based fast local level set method, Journal of Computational Physics, 155 (1999) 410-438. 28

[36] J. Gomes, O. Faugeras, Reconciling distance functions and level sets, in: Biomedical Imaging, 2002. 5th IEEE EMBS International Summer School on, IEEE, 2002, pp. 15 pp. [37] P.D. Dunning, H. Alicia Kim, A new hole insertion method for level set based structural topology optimization, International Journal for Numerical Methods in Engineering, 93 (2013) 118-134. [38] R. Courant, K. Friedrichs, H. Lewy, On the partial difference equations of mathematical physics, IBM Journal, 11 (1967) 215-234. [39] J.A. Sethian, Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science, Cambridge University press, 1999. [40] X. Guo, W. Zhang, W. Zhong, Doing topology optimization explicitly and geometrically—a new moving morphable components based framework, Journal of Applied Mechanics, 81 (2014) 081009. [41] T. Yamada, K. Izui, S. Nishiwaki, A. Takezawa, A topology optimization method based on the level set method incorporating a fictitious interface energy, Computer Methods in Applied Mechanics and Engineering, 199 (2010) 2876-2891. [42] B. Zhu, X. Zhang, S. Fatikow, Structural topology and shape optimization using a level set method with distance-suppression scheme, Computer Methods in Applied Mechanics and Engineering, 283 (2015) 1214-1239. [43] M. Otomori, T. Yamada, K. Izui, S. Nishiwaki, Level set-based topology optimisation of a compliant mechanism design using mathematical programming, Mechanical Sciences, 2 (2011) 91-98. [44] S. Wang, M.Y. Wang, Radial basis functions and level set method for structural topology optimization, International Journal for Numerical Methods in Engineering, 65 (2006) 2060-2090. [45] S. Wang, K. Lim, B.C. Khoo, M.Y. Wang, An extended level set method for shape and topology optimization, Journal of Computational Physics, 221 (2007) 395-421. [46] J. Chen, V. Shapiro, K. Suresh, I. Tsukanov, Shape optimization with topological changes and parametric control, International Journal for Numerical Methods in Engineering, 71 (2007) 313-346. [47] X. Qian, Full analytical sensitivities in NURBS based isogeometric shape optimization, Computer Methods in Applied Mechanics and Engineering, 199 (2010) 2059-2071. [48] X. Qian, Topology optimization in B-spline space, Computer Methods in Applied Mechanics and Engineering, 265 (2013) 15-35. [49] W. Zhang, W. Yang, J. Zhou, D. Li, X. Guo, Structural topology optimization through explicit boundary evolution, Journal of Applied Mechanics, 84 (2017) 011011. [50] M.Y. Wang, S. Wang, Parametric shape and topology optimization with radial basis functions, in: IUTAM Symposium on Topological Design Optimization of Structures, Machines and Materials, Springer, 2006, pp. 13-22. [51] C. Fleury, CONLIN: an efficient dual optimizer based on convex approximation concepts, Structural Optimization, 1 (1989) 81-89. [52] K. Svanberg, The method of moving asymptotes- a new method for structural optimization, International Journal for Numerical Methods in Engineering, 24 (1987) 359-373. 29

[53] K. Svanberg, MMA and GCMMA-two methods for nonlinear optimization, in, Technical Report, 2007. [54] C. Li, C. Xu, C. Gui, M.D. Fox, Distance regularized level set evolution and its application to image segmentation, Image Processing, IEEE Transactions on, 19 (2010) 3243-3254. [55] J.A. Sethian, A. Wiegmann, Structural boundary design via level set and immersed interface methods, Journal of Computational Physics, 163 (2000) 489-528. [56] M.D. Buhmann, Radial basis functions: theory and implementations, Cambridge University Press, 2003. [57] P. Wei, M.Y. Wang, Parametric structural shape and topology optimization method with radial basis functions and level-set method, in: ASME 2006 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers, 2006, pp. 463-470. [58] M.Y. Wang, The augmented Lagrangian method in structural shape and topology optimization with rbf based level set method, in: Proceedings of the Fourth China-JapanKorea Joint Symposium on Optimization of Structural and Mechanical Systems, 2006. [59] M. Mongillo, Choosing basis functions and shape parameters for radial basis function methods, SIAM Undergraduate Research Online, 4 (2011) 190-209. [60] J. Dolbow, T. Belytschko, An introduction to programming the meshless Element F reeGalerkin method, Archives of Computational Methods in Engineering, 5 (1998) 207241. [61] S. Wang, M.Y. Wang, A moving superimposed finite element method for structural topology optimization, International Journal for Numerical Methods in Engineering, 65 (2006) 1892-1922. [62] X. Xing, M.Y. Wang, B.F.Y. Lui, Parametric shape and topology optimization with moving knots radial basis functions and level set methods, 7th World Congress on Structural and Multidisciplinary Optimization, 1 (2007) 2.

30

Highlights: • A parametric level set approach for topology optimization driven by mathematical programming. • A variational method enables distance regularization inside the topology optimization loop. • Enforce the signed-distance property within a narrow band along the boundaries while keeping it flat in the rest area. • Robustness is better than conventional level set methods, especially in handling nonconvex optimization problems with multiple constraints. • The flatness of the level set function enables easier nucleation of new holes without needing topologicalm derivatives.