Topological shape optimization design of continuum structures via an effective level set method

This paper proposes a new level set method for topological shape optimization of continuum structure using radial basis function (RBF) and discrete wavelet transform (DWT). The boundary of the structure is implicitly represented as the zero level set of a higher-dimensional level set function. The interpolation of the implicit surface using RBF is introduced to decouple the spatial and temporal dependence of the level set function. In doing so, the Hamilton-Jacobi partial differential equation (PDE) that defines the motion of the level set function is transformed into an explicit parametric form, without requiring the direct solution of the complicated PDE using the finite difference method. Therefore, many more efficient gradientbased optimization algorithms can be applied to solve the optimization problem, via updating the expansion coefficients of the interpolant and then evolving the level set function and the boundary. Furthermore, the DWT is employed to handle the full matrix arising from the globally supported RBF interpolation. Several high stiffness but lightweight designs with smooth and clear structural boundaries are optimized and presented. The numerical results show that the proposed method can remarkably increase the efficiency in the topology optimization design of both the 2D and 3D structures. *Corresponding author: Mi Xiao, State Key Lab of Digital Manufacturing Equipment and Technology, Mechanical Engineering, Huazhong University of Science and Technology, 1037 Luoyu Road, Wuhan, Hubei 430074, China E-mail: xiaomi@hust.edu.cn Reviewing editor: Wei Meng, Wuhan University of Technology, China Additional information is available at the end of the article ABOUT THE AUTHORS Dr Hao Li received his PhD degree in Mechanical Engineering from Huazhong University of Science and Technology (HUST), China. He is a postdoctoral researcher with HUST since March 2016, and currently he is a visiting scholar at the University of Technology, Sydney. His research activities include structural topology optimization and design of metamaterials. Prof Liang Gao received his PhD degree in Mechanical Engineering from HUST, in 2002. Now, as a professor and head of the Department of Industrial & Manufacturing System Engineering at HUST, his research interests include operations research, optimization in design and manufacture. Dr Mi Xiao received his PhD degree in Mechanical Engineering from HUST, in 2012. He was a postdoctoral researcher in the School of Mechanical Science and Engineering (MSE) of HUST, from 2012 to 2015. Now, as a lecturer in School of MSE of HUST, his research interests include topology optimization and multidisciplinary design optimization. PUBLIC INTEREST STATEMENT An efficient level set method is proposed to solve structural topology optimization problems by using radial basis function (RBF) and discrete wavelet transform (DWT). The RBF is adopted to interpolate the discrete value of implicit level set function. Thus, the conventional level set method is converted into a relatively easier parametric form. The efficient gradient-based optimization algorithms can be directly applied to solve the optimization problem. Furthermore, the DWT is employed to compress the full interpolation matrix associated with the globally supported RBF, leading to a remarkable reduction of the CPU time under an acceptable numerical accuracy. Received: 27 June 2016 Accepted: 14 October 2016 First Published: 20 October 2016 © 2016 The Author(s). This open access article is distributed under a Creative Commons Attribution (CC-BY) 4.0 license.


PUBLIC INTEREST STATEMENT
An efficient level set method is proposed to solve structural topology optimization problems by using radial basis function (RBF) and discrete wavelet transform (DWT). The RBF is adopted to interpolate the discrete value of implicit level set function. Thus, the conventional level set method is converted into a relatively easier parametric form. The efficient gradient-based optimization algorithms can be directly applied to solve the optimization problem. Furthermore, the DWT is employed to compress the full interpolation matrix associated with the globally supported RBF, leading to a remarkable reduction of the CPU time under an acceptable numerical accuracy.

Introduction
Topology optimization of continuum structures has been recognized as one of the most challengeable design problems in the field of structural optimization. Topology optimization is essentially a numerical process to iteratively re-distribute a given amount of materials inside the design domain to determine the best material layout of the structure subject to appropriate boundary conditions, so as to improve structural performance until the objective function is optimized under design constraints. Over the past two decades, the topology optimization (Bendsøe & Sigmund, 2003) has been applied to a range of engineering fields, including structures, mechanisms and materials. Several methods have been developed for topology optimization problems.
After one of the pioneers' work in 1988 (Bendsøe & Kikuchi, 1988), the topology optimization of continuum structures has experienced considerable development by using different numerical analysis approaches, such as the finite element method (FEM). In the work of Bendsøe and Kikuchi (1988), the microstructures with regularly shaped holes are positioned in the design space, and the topology optimization problem is transformed into a size optimization problem of the microstructures by using the optimality criteria (OC) method (Zhou & Rozvany, 1991). To improve the efficiency of the topology optimization, a generalized homogenization method, the solid isotropic microstructure with penalization (SIMP) approach (Rozvany, Bendso̸ e, & Kirsch, 1995), has been developed. Other than the homogenization method which employs the size, position and angle of a microstructure as the design variables, the SIMP scheme introduces the artificial densities of the elements to denote the material properties of the elements. The SIMP method has gained its popularity due to the conceptual simplicity and implementing easiness. Both the SIMP and homogenization methods are mainly based on the FEM, to relax the original discrete topology optimization problem as the continuous problem, which allows the occurrence of the intermediate values of the design variables. Another method is known as evolutionary structural optimization (ESO) (Xie & Steven,1993), which usually treats the topological design of the structure as a combinational optimization problem with discrete design variables, and removes the invalid material from the domain iteratively based on pre-determined criteria. A family of pointwise density-based interpolation (PDI) methods (Kang & Wang, 2011;Luo, Zhang, Wang, & Gao, 2012) is also developed for topology optimization of structures. However, most of the above topology optimization methods are commonly characterized with zigzag geometrical boundary and blur material interface, which makes it difficult to extract the geometry of the final design.
Since Sethian and Wiegmann (2000) initially introduced the standard level set method into the area of structural optimization, this method has been recognized as an alternative approach for the topological shape optimization design of structures with smooth and clear boundaries (Allaire, Jouve, & Toader, 2004;Wang, Wang, & Guo, 2003). Despite its favorable merits (Allaire et al., 2004;Wang et al., 2003), the standard level set method has several numerical difficulties in handling the Hamilton-Jacobi PDEs, e.g. the re-initialization of the level-set function, the extension of the velocity field and the restriction of the CFL condition (Allaire et al., 2004;Luo, Wang, Wang, & Wei, 2008;Wang et al., 2003). To avoid the unfavorable numerical features, the parametric level set methods (Luo et al., 2008;Wang, Luo, Zhang, & Kang, 2014;Wang & Wang, 2006) have recently been developed for a variety of topology optimization problems.
Radial basis function (RBF) can be efficiently utilized in interpolating scattered data within the higher-dimensional spaces (Wendland, 2005). This capacity makes them very attractive in the level set methods. In the work of Wang and Wang (2006), the multiquadric and inverse multiquadric splines were utilized as the globally supported RBF (GS-RBF) to interpolate the level set function. The implicit level set model was parameterized to avoid the direct solution of the Hamilton-Jacobi PDEs. The positive features of GS-RBFs include the solvability, global smoothness, accuracy and convergence of the interpolation problem. However, further application to topology optimization problems of continuum structures was restricted by the prohibitive computational cost due to the global support. Moreover, the level set front propagation was obtained by solving a nonlinear set of ordinary differential equations (ODEs) using the steepest descent algorithm, rather than more efficient gradient-based optimization algorithms. To overcome the shortcomings of GS-RBF interpolants, Luo et al. (2008) incorporated the compactly supported RBF (CS-RBF) into the framework of parametric level set method, which is more efficient and effective for topology optimization problems, although the numerical accuracy is lower than GS-RBF. Another parametric level set method (Ho, Wang, & Zhou, 2012) was developed by using a dynamic moving RBF knots scheme and the Partition of Unity (POU) method to alleviate the high computation cost due to the fully dense matrix in the parameterization using GS-RBF. This paper will propose another new parametric level set method to rich the family of alternative level set methods for topology optimization of structures. It is well known that evaluating the linear system arising from RBFs normally requires O(N 3 ) flops and O(N 2 ) memory, which means the computational cost is considerably expensive and the optimization may be failed when handling the largescale structural optimization problems. Indeed, a sparse system is a reasonable solution to deal with this numerical difficulty. Therefore, in this paper, we will propose a more effective and efficient level set method to overcome the shortcomings while inheriting the advantages of the existing parametric level set method. The GS-RBFs will be employed to interpolate the level set function, while DWT will be incorporated into the framework of the parameterization level set method, to provide a sparse linear system for the interpolation of discrete level set function. The DWT is a multi-resolution decomposition for input data and has been widely used in signal processing, image compression, computer vision, denoising, face detection and so on (Hsia, Guo, & Chiang, 2012). Apart from these successful applications, it has recently been used for handling the dense and fully populated matrices (Chen, 1999;Ravnik, Škerget, & Hriberšek, 2004), because it is capable of rapidly capturing the critical data with only a few useful coefficients among scattered information set. In another word, this strong capability is mainly based on the fact that the majority of data are correlated with time and frequency. It is expected that the percentage of zero elements in the matrix can be over 90%.
Several numerical examples are used to demonstrate the effectiveness of the proposed level-set method.

Implicit boundary representation
Regarding the level set method, the design boundary of the structure is often embedded implicitly into a higher-dimensional scalar function. In doing this, the merging and breaking of the moving boundary during iterations can be easily described or tracked. Let us define D ⊂ ℝ d (d = 2, 3), which is the reference domain containing all admissible shapes of Ω. According to the level set method, each part in the design domain can be described as (Sethian & Wiegmann, 2000): Introducing the pseudo time t to enable the dynamic motion of the level set function, differentiating both sides of Φ(x, t) respect to t and applying the chain rule, the propagation of the level set function can be mathematically defined as the following Hamilton-Jacobi equation: where V n is the normal velocity that can be expressed as: Now, the key of the topological shape optimization of the structure is to find the feasible solution of the Hamilton-Jacobi PDE. However, as mentioned previously, the direct solving of Hamilton-Jacobi PDE is not an easy task and will hinder the application of the standard level set method associated with the topological shape optimization problems.

Parametric model based on RBF
In order to overcome the unfavorable numerical issues in the conventional level set methods, the GS-RBF will be used to implement the interpolation of the level set function. In these formulations, the level set function can be approximated by centrally positioning a series of radially symmetric RBFs at the knots over the design domain. Thus, the original time-dependent level set function Φ(x, t) is decoupled by a group of expansion coefficients and RBF interpolants. The level set-based topological shape optimization problem is then transformed into a parametric optimization problem, i.e. a size optimization of the expansion coefficients, which can be solved by the well-established gradient-based approaches in the space of the parameters.
With the RBF interpolation, the level set function can be approximated as follows: where the vector of shape functions are defined by and the vector of expansion coefficients are given as: The Gaussian RBFs (Wendland, 2005) are used due to its simplicity, positive definiteness and high performance: where r is calculated by the Euclidean norm on ℝ d : and where x i denotes the position of interpolant knot.
The Gaussian RBFs are among the family of GS-RBFs, which have many attractive features, such as high level of interpolation accuracy and smoothness. However, two crucial issues may prevent the further application of these GS-RBFs from topological shape optimization of the structure. The first is how to select the free shape parameter c in (7), which has an obvious influence upon the accuracy of interpolation (Wang & Wang, 2006). The second is the computational efficiency as well as the computer storage caused by the fully dense matrix, especially for large-scale optimization problems (Luo et al., 2008). The latter will be studied and illustrated in the following sections, while the former one has already been discussed in studies (Wendland, 2005) and it is beyond the scope of this paper.
Substituting (4) into (2), the Hamilton-Jacobi PDE will change to a parametric form with the timespace independent variables as follows: Thus the normal velocity field is given as follows: It can be seen that the normal velocity is naturally extended to all the RBF knots over the design domain, since φ(x), α(t) and ̇ (t) are all evaluated over the knots in the design space. To facilitate the discussion, we introduce the following theoretically invertible matrix A: Substituting (11) into (4), the level set values over the knots can be rewritten in the form of the following linear equations: Since the level set values Φ 0 are pre-known in the initialization phase of the optimization, the initial values of the expansion coefficients can be obtained by solving the linear system in (12): The interpolation is automatically applied to the knots over the entire design domain rather than only the points on the boundary. For the numerical simplicity, the identical grids are used to implement FEA calculation and RBF process.

Matrix compression using DWT
In this paper, the DWT (Chen, 1999;Ravnik et al., 2004) is implemented to improve the efficiency of the numerical process for the RBF-based parametric model. Wavelet transform is a mathematical tool developed especially for saving the computational time and memory requirement. By utilizing the invertible transformation and an appropriate threshold, it permits sparse coefficients matrices for interpolation at a low computational cost, and the approximate solution is obtained by inversely transforming the resultant vector using the wavelet basis into a new one in the original basis. More importantly, this matrix compression technique can act as a black box to be incorporated into the established parametric formulation. This means it can provide a quite sparse linear system through the DWT strategy without complicated modification to the existing model, and only very few amount of extra CPU time will cost by specified operators of the transform.
A pyramidal scheme is introduced to transform the original basis into a wavelet basis (Beylkin, Coifman, & Rokhlin, 1991). We adopt the Haar wavelets that are the compact support wavelets of Daubechies with one vanishing moment, since that they have relatively simple conception, constant scaling function and non-overlapping support (Daubechies, 1988). Regarding the linear system given in (12), the vector or can be transformed into the wavelet basis recursively as: where j = 1, 2,…,n are levels of the decomposition, and k = 1, 2,…,n = 2 j denote the numbers of components in the vector. s j k and d j k are elements in sub-vectors and can be regarded as periodic sequences with the period 2 n−j . It is remarked that s 0 k represent the components in the original vector. For the Haar wavelet, we have M = 1, which means this wavelet is characterized by 2 wavelet filter coefficients (h i , g i ) (Ravnik et al., 2004).
The decomposition of matrix A is actually a 2D wavelet transform for matrix, which is composed of two 1D DWTs for vectors. A typical one level and 2D wavelet transformation can be summarized as in Figure 1. The matrix is initially divided into two sub-matrixes with the same size. The results are then preceded with two vertical 1-D DWTs separately. An operation called downsampling is applied to the filtered results. The pair-wise filters, involving a low-pass filter and a high-pass filter, eventually decompose the matrix into the low-low (LL), low-high (LH), high-low (HL) and high-high (HH) wavelet frequency bands with a quarter size of the initial matrix, respectively. The requirement that the number of components in a vector must be a power of 2 when conducting the fast wavelet transform in the traditional DWT application has been extended to an arbitrary vector length (Ravnik et al., 2004).
For a transformed vector, the re-construction process to inversely convert it into a vector with the original basis is derived as follows: For the simplicity, the orthogonal matrix W (Chen, 1999) is introduced to represent the DWT process given in Equations (14-17), alternatively. Regarding the linear system expressed in (12), the ̄= and ̄ = are the wavelet forms for the vectors Φ and α. Similarly, the matrix A in the original basis is then represented by ̄ = −1 as a relatively new matrix in the wavelet basis. Only a few significant coefficients in ̄ are essential for manipulating the entire data-set after the DWT. This means we can compress ̄ to build a sparse matrix ̃ by sweeping out the noise elements and only conserving the entries which are greater than an appropriate tolerance threshold (Chen, 1999;Ravnik et al., 2004). With the pre-multiplication of W on both sides of Equation (12), the linear system can be rewritten as: Substituting the transformed vectors ̄ , ̄ as well as the compressed matrix ̃ into (18), we have a new sparse linear system: This system is easy to be handled. Solving the equations to obtain ̄ or multiplying two components to compute the product ̄ will require relatively lower computational cost. Eventually, the expansion coefficients are calculated by the inverse DWT as = −1̄ , and the level set function values are obtained by = −1̄ . By far, the DWT-based matrix compression technique has been incorporated into the RBF-based parametric method, which will lead to a linear system with lower computational cost as well as computer storage.

Design sensitivity analysis
For linear elastic structural problem, the topology optimization problem of stiffness design under a global volume constraint is formulated with the level set-based implicit representation: The structural compliance is f (u, u) = T ij (u)C ijkl kl (u)∕2. V max is the prescribed volume constraint, u is the displacement field and v is the virtual displacement field. u 0 is defined on the admissible Dirichlet boundary ∂ Ω. H(Φ) is the Heaviside function.
The energy bilinear functional a (u, v, Φ) and the load linear form l (u, v) in the state equation can be written in their weak form, respectively, as follows: where ɛ is the strain field, p is the body force and τ is the boundary traction. δ(Φ) represents the partial derivative of the Heaviside function, i.e. Dirac function (Allaire et al., 2004;Wang et al., 2003).
To solve the constrained optimization problem, the Lagrange multipliers λ and Λ are introduced to transform it into an unconstrained one, as follows: In our approach, we use shape derivative analysis to achieve the design sensitivities. According to Wang et al. (2003), the shape derivative of the Lagrange function can be obtained by: where β(u, u, Φ) is the so-called shape gradient density, readers can refer to the relevant literatures (Allaire et al., 2004;Luo et al., 2008;Wang et al., 2003) for more details.
Substituting normal velocity defined in (10) into (24), the shape derivative of the Lagrange function can be rewritten by Alternatively, the shape derivative can be also derived by the chain rule: Therefore, by comparing the corresponding terms in (25) and (26), design sensitivities for the objective function and constraint can be determined as: Here the volume integration is adopted to improve the computational efficiency in numerical process, due to the fact that it doesn't need to calculate along the design boundary (Wang et al., 2014).

Optimization algorithm
As discussed previously, the parametric level set formulation can be solved by many well-established gradient-based methods. Here, we adopt the OC method to update the design variables, because of its implement simplicity and high efficiency in dealing with large-scale optimization problems with only a few constraints (Rozvany et al., 1995;Zhou & Rozvany, 1991). The OC-based method implements an explicit heuristic scheme based on the Kuhn-Tucker conditions for the constrained optimization problem to update the design variables. According to the work by Luo et al. (2008), the updating scheme for

Numerical examples
In this section, several examples are used to showcase the characteristics of the proposed method with different settings. For the sake of simplicity, the coordinates of the RBF knots is supposed to be identical with the finite element nodes, and the shape parameter of the Gaussian kernel is equivalent to the size of the finite element. The moving limit is 0.003, and the damping factor is 0.3 for the OC method, respectively. The optimization algorithm is terminated when the relative difference of the objective function values within the latest two iterations is less than 0.0001. All the routines are done with MATLAB and operated on a computer with a CPU 2.67 GHz processor as well as a 4GB RAM.

2D cantilever beam
As shown in Figure 2, the design domain is a 2D cantilever beam with the length and width ratio of 2:1. A concentrated force F = 1 is applied at the center of the right edge, and the left side of the beam is fixed. Regarding the artificial material model, the Young's modulus E = 1 for solid materials, E = 10 −3 for weak materials and Possion's ratio v = 0.3. The aim is to minimize the mean compliance of the structure for a stiffness design, and the maximum allowable material usage inside the design domain is 50%.
In this example, we first study the effect of the threshold in the proposed method. The cantilever beam is discretized with 100 × 50 = 5,000 rectilinear finite elements. Accordingly, the interpolant matrix ̄ will have a large number of coefficients (101 × 51) 2 . In the proposed approach, one of the keys is to threshold noise elements in the system matrix ̄ , which is constructed by using the wavelet basis. Here, we employ the thresholding operation proposed by Ravnik et al. (2004): where ā ij denotes the single element in ̄ , and m is calculated by the average value of absolute coefficients in this system matrix. It is convenient to change the thresholding limit with factor κ ≥ 0.
As shown in Figure 3, the shape and topology of the boundary evolve toward the optimality with κ = 1. Figure 4 displays the changes for the corresponding level set surface. It can be found that the proposed approach can maintain the advantages of handling shape fidelity and topological changes in the level set-based methods. The optimal design in Figure 3(f) has a smooth boundary, similar to the results given in the literatures (Allaire et al., 2004;Luo et al., 2008;Wang et al., 2003). Convergent histories of the objective function and volume constraint are depicted in Figure 5. It takes 309 steps to minimize the mean compliance to 60.5706, in which the first 5 iterations are mainly utilized to push back the volume constraint due to the violation. It is noticed that the convergence is stable and the volume constraint is conservative, which demonstrate the effectiveness of the proposed method.
0, otherwise.   We further examine the optimization results with different factor κ. As a comparison, the optimization problem is solved by the original parametric method without including the DWT and the proposed method with the DWT, respectively. For each case, the design problem and the relative optimization parameter settings are identical. We only apply the different κ values in the DWT. Figure 6 displays the optimal results of different cases. It is obvious that the topologies for all the designs are nearly the same. Relative statistics for the tested cases are illustrated in Table 1. Notice that the main computational cost include FEA and optimization process (involve sensitivity analysis). FEA time is fixed if we don't change the meshes. Thus, we only need to compare the optimization times to showcase the efficiency of the algorithm. In addition, to start the iteration, one must first solve (19) to initialize the expansion coefficients, which is called as "fitting" process and will influence the CPU time as well.
The evidence implies that with the larger κ values, the sparser system will be observed. However, this does not apply any remarkable impact on the optimal compliances. Obviously, we see that the computational cost and the computer storage are dramatically reduced via the implementation of the proposed method. Nonetheless, It should be pointed out that the fitting times for Cases 1(b) and 1(c) are larger than that of Case 1(a), due to that superabundant processing time are cost to complete the wavelet transform when the interpolant matrices are not enough sparse.

3D Supported structure
The supported structure of a satellite with a dimension of 1 m × 4.8 m × 4.8 m is given in Figure 7. Regarding the whole design domain, four vertical forces F = 500 kN are applied at the top surface and nine points are fixed on the bottom face. The Young's modulus E = 44.8 GPa for solid materials, the Possion's ratio v = 0.3, and the material density ρ = 1.78 g/cm 3 . The aim is to minimize the strain energy, and the maximum allowable material usage is 15%. The design domain is meshed with 10 × 48 × 48 eight nodes finite elements, and only 10 × 24 × 24 elements are taken into considered because of the symmetry. The proposed parametric formulation is used, with a thresholding factor κ = 1.
With the proposed level-set method, the percentage of zero elements in the interpolation matrix A is 99.28%. The iterative process is shown in Figure 8 with the minimum strain energy 278.0171 after 159 steps. After topological changes and shape fidelity, the optimal design is achieved with a smooth and distinct boundary, which can be conveniently fabricated by the additive manufacturing technology. The optimal result in Figure 8(f) is similar to the result given in Pan, Chen, and Wang (2004), which has theoretically proven to be the optimal structure for supporting. Figure 9 shows the convergence of the objective function and the volume constraint. It can be seen that the violated volume constraint in the initial design stage is pushed back after the first 15 iterations, and is conserved in the following iterations. Furthermore, we remark that the proposed method can give a stable optimization process with acceptable CPU time (i.e. 79.9375 s per step for the optimization and 2.0156 s for the fitting), as well as the memory requirement.

Conclusions
This paper has proposed an effective topology optimization method for continuum structures by combining the favorable features of RBF and DWT. The RBF is used to parameterize the implicit level set surface, and uncouple the Hamilton-Jacobi PDE into a set of algebraic equations. Consequently, the complex topology optimization problem is converted into a size optimization problem of the expansion coefficients of the interpolant. In the level set formulation, the X-FEM scheme provides a more accurate way to calculate the strains when the level-set boundary crossing finite elements. The fully dense interpolant matrix of the GS-RBFs often leads to highly cost CPU time and expensive memory requirements. Therefore, the DWT approach can compress the matrix and produce a sparse interpolation system, which can greatly improve the computational efficiency. The numerical results show that the thresholding parameter κ plays a significant role in the proposed method. A smaller value will reduce the efficiency of the optimization process, while a larger value will cause the loss in interpolation accuracy and iterative stability. Hence, numerical experience more or less will be