1 Engineering Faculty, China University of Geosciences, Wuhan 430074, China
2 Hubei Engineering University, Xiaogan 432000, China
To cure imperfections such as low accuracy and the lack of ability to nucleate hole in the conventional level set-based topology optimization method, a novel method using a trapezoidal method with discrete design variables is proposed. The proposed method can simultaneously accomplish topology and shape optimization. The finite element method is employed to obtain element properties and provide data for calculating design and topological sensitivities. With the aim of performing the finite element method on a non-conforming mesh, a relation between the level set function and the element densities field has to be clearly defined. The element densities field is obtained by averaging the Heaviside function values. The Lagrange multiplier method is exploited to fulfill the volume constraint. Based on topological and design sensitivity and the trapezoidal method, the Hamilton-Jacobi partial differential equation is updated recursively to find the optimal layout. In order to stabilize the iterations and improve the efficiency of the algorithm, re-initiation of the level set function is necessary. Then, the detailed process of a cantilever design is illustrated. To demonstrate the applications of the proposed method in bridge construction, two numerical examples of a pylon bridge design are introduced. It is shown that the results match practical designs very well, and the proposed method is a helpful tool in bridge design.
Keywords: Level set method, Sensitivity analysis, Topology optimization, Trapezoidal method.
open-access license: This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International Public License (CC-BY 4.0), a copy of which is available at: https://creativecommons.org/licenses/by/4.0/legalcode. This license permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
* Address correspondence to this author at the Engineering Faculty, China University of Geosciences, Wuhan 430074, China; Tel: 86-02767883428; Fax: 86-02767883507; E-mail: email@example.com
Topology optimization is conducive to the conceptual design of a structure and enables engineers to find optimal material distribution. A wide range of fascinating applications of topology optimization has been developed [1J.N. Richardson, S. Adriaenssens, P. Bouillard, and R.F. Coelho, "Multiobjective topology optimization of truss structures with kinematic stability repair", Struct. Multidiscipl. Optim., vol. 46, no. 4, pp. 513-532, 2012. [http://dx.doi.org/10.1007/s00158-012-0777-5] -3K. Maute, and G.W. Reich, "Integrated multidisciplinary topology optimization approach to adaptive wing design", J. Aircr., vol. 43, no. 1, pp. 253-263, 2006. [http://dx.doi.org/10.2514/1.12802] ]. Sophisticated methods have emerged to help engineers deal with the topology optimization problem of discrete structures. However, there are more complexities by far in the topology optimization problem concerning a continuum structure. The methods available for engineers to solve continuum optimization problems are still limited and imperfect. Accordingly, research on the continuum optimization problem is still undertaken for better optimization results and accuracy.
After almost 30 years of development, four crucial methods have stood out for the topology optimization of the continuum structure, namely, the homogenization method [4A.R. Díaz, and M.P. Bendsøe, "Shape optimization of structures for multiple loading conditions using a homogenization method", Struct. Optim., vol. 4, no. 1, pp. 17-22, 1992. [http://dx.doi.org/10.1007/BF01894077] ], the level set-based methods [5Z. Luo, N. Zhang, J.C. Ji, and T. Wu, "A meshfree level-set method for topological shape optimization of compliant multiphysics actuators", Comput. Methods Appl. Math., vol. 223–224, pp. 133-152, 2012.], the evolutionary structural optimization (ESO) approach [6X. Huang, and Y.M. Xie, Evolutionary Topology Optimization of Continuum Structures: Methods and Applications., John Wiley & Sons: New York, 2010. [http://dx.doi.org/10.1002/9780470689486] ], and the solid isotropic material with penalization (SIMP) scheme [7M.P. Bendsøe, "Optimal shape design as a material distribution problem", Struct. Optim., vol. 1, no. 4, pp. 193-202, 1989. [http://dx.doi.org/10.1007/BF01650949] ].
The homogenization method originated from the pioneering work of Bendsoe and Kikuchi in the late 1980s [8M.P. Bendsoe, and N. Kikuchi, "Generating optimal topologies in structural design using a homogenization method", Comput. Methods Appl. Mech. Eng., vol. 71, pp. 197-224, 1988. [http://dx.doi.org/10.1016/0045-7825(88)90086-2] ]. In the 1990s, most of the topology optimization work depended on the homogenization method [9G. Allaire, Shape Optimization by the Homogenization Method., Springer Verlag: New York, 2001., 10G. Allaire, E. Bonnetier, G. Francfort, and F. Jouve, "Shape optimization by the homogenization method", Numer. Math., vol. 76, pp. 27-68, 1997. [http://dx.doi.org/10.1007/s002110050253] ]. In addition, the SIMP scheme can be regarded as a variant of the homogenization method [11M.P. Bendsøe, and O. Sigmund, "Material interpolation schemes in topology optimization", Arch. Appl. Mech., vol. 69, pp. 635-654, 1999. [http://dx.doi.org/10.1007/s004190050248] , 12M.P. Bendsøe, and O. Sigmund, Topology Optimization: Theory, Methods, and Applications., Springer-Verlag: Berlin, Heidelberg, 2004. [http://dx.doi.org/10.1007/978-3-662-05086-6] ]. In the late 1990s and early 2000s, the SIMP scheme supplanted the homogenization method and became popular for researchers due to its simplicity. The basic idea of ESO was developed from the one-way hard-kill strategies. The structure is first discretized, and the material of elements that have the lowest strain energy density are set to be null in each iteration [13O. Sigmund, and K. Maute, "Topology optimization approaches, a comparative review", Struct. Multidiscip. O., vol. 48, no. 6, 2013.].
Moreover, researchers ingeniously innovated some multi-discipline methods [14M. Ohsaki, "Genetic algorithm for topology optimization of trusses", Comput. Struc., vol. 57, no. 2, pp. 219-225, 1995. [http://dx.doi.org/10.1016/0045-7949(94)00617-C] , 15K. Suresh, "Efficient generation of large-scale pareto-optimal topologies", Struct. Multidiscipl. Optim., vol. 47, no. 1, pp. 49-61, 2013. [http://dx.doi.org/10.1007/s00158-012-0807-3] ]. The level set method was first created by Sethian and Osher [16S. Osher, and J.A. Sethian, "Front propagating with curvature-dependent speed: Algorithms based on Hamilton–Jacobi formulations", J. Comput. Phys., vol. 79, no. 1, pp. 12-49, 1988. [http://dx.doi.org/10.1016/0021-9991(88)90002-2] ]. In the early 2000s, the level set-based method was introduced in the topology optimization discipline and became increasingly popular. Level set-based topology optimization methods could simultaneously evade the checkboard problem and accomplish topology and shape optimization, which are unreachable by density-based methods. The boundary is characterized by the principal process of the level set method with the Hamilton-Jacobi equation, making the boundary evolve and allowing boundary manipulation. The evolution of the boundary uses the analytical design sensitivity. Therefore, the design sensitivity always corresponds to the boundary velocity.
Topology changes of structures can be made easily by the level set method. However, they will encounter trouble when inadequate holes are initially laid within the structure for their inability to increase the number of holes, which is a direct consequence of the maximum principle. Hence, the speculation regarding the number of holes is of great importance for the conventional level set methods for topology optimization. To overcome the problem, the topological sensitivities could be introduced when the level set function is changed, facilitating the nucleation of holes [17G. Allaire, F.D. Gournay, F. Jouve, and A.M. Toader, "Structural optimization using topological and shape sensitivity via a level set method", Control Cybern., vol. 34, no. 1, pp. 59-80, 2005.]. The topological sensitivities are acknowledged to be one of the most important tools in topology optimization.
Using the topological sensitivities, this study proposes a more accurate level set-based topology optimization method. In place of the tapping upwind scheme, a method merely assures first order accuracy, the trapezoidal method [18G. Strang, Computational Science and Engineering., Wellesley-Cambridge Press: Massachusetts, 2007.] can be utilized to promote boundary changes of the structure, bringing about more accurate results. The topological sensitivities and the design sensitivities are combined to establish the velocity of boundary. Then, the trapezoidal method is used to solve the Hamilton-Jacobi equation.
2. LEVEL SET METHOD
2.1. Level Set Function
A free moving boundary that encircles the design space is set to be the zero level set. This is the underlying principle of the level set method. Any change occurring in the boundary is susceptible to the evolvement of the level set function, which occupies a higher order one-dimensional space [19S. Osher, and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces., Springer-Verlag: New York, 2003. [http://dx.doi.org/10.1007/b98879] -21M. Burger, "A framework for the construction of level set methods for shape optimization and reconstruction", Interface Free Bound, vol. 5, no. 3, pp. 301-330, 2003. [http://dx.doi.org/10.4171/IFB/81] ]. As Fig. (1) shows, D denotes the structural design domain, which is the subject of the topology optimization. The material domain is denoted by Ω, which is specified using the boundary
Fig. (1) Design domain and material domain.
Fig. (2) illustrates that the two-dimensional boundary is defined by a zero level set of signed distance function, which occupies a three-dimensional space and is also a level set function. When the structure is a subset of some domain Ω, which lies in a bigger design domain D, the subsequent properties manifest themselves as follows:
Level set function and zero level set.
, as before, indicates the lower dimensional boundary, while Ф(x) denotes the level set function. In the optimization of structure, for simplicity, the level set of zero is usually used to indicate the boundary as
Operating the preceding formula with differentiation gives the evolution formula for the level set function, which is defined with the subsequent Hamilton-Jacobi equation:
where denotes normal component of velocity, while t denotes pseudo time.
When a function, which is described with less dimensions, is embedded into a function, which is described with more dimensions, there will be many favorable features [22N.P. Van Dijk, G.H. Yoon, F. Van Keulen, and M. Langelaar, "A level-set based topology optimization using the element connectivity parameterization method", Struct. Multidiscipl. Optim., vol. 42, no. 2, pp. 269-282, 2010. [http://dx.doi.org/10.1007/s00158-010-0485-y] ]. Furthermore, representing the function described with more dimensions by a signed distance function facilitates the update process. Hence, the signed distance function is frequently used to represent the level set function at first. However, the degeneration of convergence properties emerges in the process of updating. Hence, the re-initialization techniques are required to regulate the level set function [23N.P. Van Dijk, M. Langelaar, and F. Van Keulen, "Explicit level-set-based topology optimization using an exact heaviside function and consistent sensitivity analysis", Int. J. Numer. Methods Eng., vol. 91, no. 4, pp. 67-79, 2012. [http://dx.doi.org/10.1002/nme.4258] ].
2.2. Structure Discretization
Prior to the optimization process, the structure is required to be discretized, after which the analytical sensitivities are calculated. With the aim of performing the finite element method on a non-conforming mesh, the relation between the level set function and the element density field pe, which is usually obtained by averaging the Heaviside function H, must be clearly defined:
in which Ωe is the material domain of e.
An approximate Heaviside function could be utilized to evaluate the element density, for example, a third order polynomial:
where ε is the lower bound for the element densities and h denotes the bandwidth of the approximate Heaviside. When h → 0, it is transformed to the exact Heaviside and physically equals the volume fraction of material in each cell.
The density evaluated by the equation (5) is eligible to be directly applied to assess the stiffness, which is called the ersatz material approach. The stiffness element is evaluated by the following equation:
in which Ke denotes the elemental initial stiffness matrix, and u denotes displacement vector, while f denotes the external loading vector.
3. FORMULATION OF TOPOLOGY OPTIMIZATION PROBLEM
The objective of optimization is to reduce the weight of the structure and simultaneously find the optimal distribution of materials to maximize the structural stiffness. Therefore, structural compliance is an objective function. h(pe, u) Linear constraint is employed to demonstrate the linear attribute that the structure demonstrates, and nonlinear constraint g(pe) is used to demonstrate the volume fraction. The mathematical model of the topology optimization problem for continuum is defined by:
where f is the external loading vector, Vmax and is the desired volume fraction. The element density
varies between the lower bound and ε 1. The definition of ε is to evade singularity. The element density
is evaluated by equation (5).
4. SENSITIVITY ANALYSIS
4.1. Design Sensitivity
The evolvement of the level set function demands the message from the element sensitivities. If the moving boundary satisfies the traction-free boundary conditions, the subsequent formula could be used to define the design sensitivity of the compliance objective
is a coefficient, and it is a result of equation (7) for the avoidance of singularity, while u denotes the global displacement vector, and ue denotes the element displacement vector [24F. Van Keulen, R.T. Haftka, and N.H. Kim, "Review of options for structural design sensitivity analysis. Part 1: Linear systems", Comput. Methods Appl. Math., vol. 194, no. 30–33, pp. 3213-3243, 2005.]. There is a significant difference between the formula of design sensitivity and the one employed by SIMP, as the densities defined in this work have been polarized while the avoidance of intermediate densities penalization is utilized. Straightforwardly, the design sensitivity with respect to volume constraint g(pe) is identical to that in [25V.J. Challis, "Multi-Property Topology Optimisation with the Level-Set Method,", Ph.D. thesis, The University of Queensland, Brisbane St Lucia, QLD, Australia, 2009.]:
4.2. Topological Sensitivity
Let Ω be an open-bounded domain with a smooth boundary and set Ωe to be a new created domain, then:
in which the boundary can be expressed by:
where denotes a sphere of radius ε with its center located at
[26A.A. Novotny, and J. Sokołowski, Topological Derivatives in Shape Optimization., Springer-Verlag: Berlin, Heidelberg, 2013. [http://dx.doi.org/10.1007/978-3-642-35245-4] ]. Hence, the new domain Ωe, which has the tiny hole
, emerges. Usually a shape functional Ψ(Ωe), which has association with the topologically perturbed domain, accepts subsequent topological asymptotic expansion, which is expressed by:
in which Ψ(Ω) denotes a shape functional belonging to the unperturbed domain. Moreover, f(ε) denotes a function whose values are always positive, so that f(ε) → 0 when ε → 0. Then, function
, which could be considered a first order correction of Ψ(Ω) to approximate, Ψ(Ωε) is used to express the topological sensitivity of ψ at
. Rearranging equation (17) gives the following:
Limit passage ε → 0 of the preceding equation results in:
An inference about the topological optimization could be made by exploiting its geometrical interpretation. However, for the sake of making good use of topological sensitivity, an explicit numerical formula should be given. Topological sensitivity with respect to the two-dimensional objective function
, which aims to minimize the structural compliance and satisfies the traction-free boundary conditions on the newly created hole as regards a sphere with a radius of one as the model hole, is
where ueT(KTr)e ue is the finite element approximation to the product tr(σ)tr(ε*) in which σ and ε* separately denote the stress tensor and strain tensor [27G. Allaire, F. Jouve, and A.M. Toader, "Structural optimization using sensitivity analysis and a level-set method", J. Comput. Phys., vol. 194, no. 1, pp. 363-393, 2004. [http://dx.doi.org/10.1016/j.jcp.2003.09.032] ]. In (17) λ and µ denote the Lamé parameters. The topological sensitivity with respect to volume constraint
, which involves a sphere with a radius of one, is expressed by:
5. PROPOSED METHOD
For the compliance minimum problem, we propose an iterative algorithm incorporated with the level set method, design sensitivity, and topological sensitivity. We have implemented the algorithm with MATLAB. A level set function is established according to material boundary. As discussed, it is preferable to choose the signed distance function as a level set function for its virtue of simplicity shown in the function updating process. Then, the state equations are discretized on finite elements. For simplicity, bi-linear, four-node, square elements are used. After the process of finite element analysis, the global displacement vector u and elemental displacement vector ue will be returned, and the algorithm updates the level set function.
In the method proposed in this work, the sign of the design sensitivity is reversed and employed as the normal velocity vn for the boundary of the design. For the sake of avoidance of a local minimum, which consists of insufficient holes, topological sensitivity is incorporated into the proposed update scheme. A revised trapezoidal method is utilized to solve formula (3) of the update scheme so that a more accurate result is obtained. In comparison with the conventional level set method, the forward Euler method is always exploited. Based on formula (3) and the trapezoidal method, a revised and novel update scheme for the level set function is defined as follows:
in which is a positive parameter, and Δt is a step size of time. To strike a balance between the influence of design sensitivity
and that of topological sensitivity must be well defined.
In addition, vn and η are the velocity terms that are defined by design sensitivity and topological sensitivity. To meet the volume constraint, both terms are defined by the Lagrangian design sensitivity and topological sensitivity:
are parameters. They are updated at each iteration and described through the subsequent formula:
where Vk is the volume of the current optimized result after k iterations. The coefficient
should be a static coefficient. With Lagrangian and sensitivity outcomes, velocity vn, which is defined within element at k iterations, is expressed by:
The purpose of the velocity term η, which is chosen based on topological sensitivity, is to nucleate holes. However, nucleating a hole in a void region is pointless. Hence, g is defined by:
in which δTL denotes the Lagrangian topological sensitivity that is expressed using:
Then, an iteration of optimization has finished. If the results converge, the whole optimization process will terminate. Otherwise, another iteration will begin.
The property of the level set function deteriorates as the update of the level set function proceeds. Thus, re-initialization of the set level set function back to the signed distance function is critical and required after a couple of steps. Without re-initialization, the optimized result would become unacceptable.
Finite element analysis is an extremely time-consuming process. Usually, the changes in the result of finite element analysis are trivial within a couple of update steps. With this observation, a parameter is set to decide how often finite element analysis is performed, and it results in a substantial decrease in the whole time of optimization. The flow chart of the proposed method is shown in Fig. (3).
Flow chart of proposed method.
6. NUMERICAL EXAMPLES
6.1. Design of a Cantilever
We consider the design of a cantilever that has a span of 3000 mm, height of 1000 mm, Young’s modulus = 100 GPa, and Poisson’s ratio v = 0.3. To simplify the problem, the material is deemed to be isotropic homogeneous. As depicted in Fig. (4), the design domain is subject to a loading of 100 N at the right center and is fixed on the left side. The volume fraction is 50%.
Design domain of a cantilever.
The whole design domain is discretized by a 60×20 finite element mesh. The original design domain, midway design results, final optimal results, and their corresponding level set function are shown in Fig. (5). Initially, the zero-contour of the level set function is set to be the whole design domain/material domain. The internal holes, which have to be set initially in conventional level set topology optimization, are not embedded. Then, the material domain shrinks due to the influence of the volume constraint. At the same time, the algorithm endeavors to nucleate internal holes spontaneously without initial assistance. After a couple of steps, the structure reaches its optimal result.
Zero level set contours and level set surfaces (60×20) of a cantilever.
6.2. Design of Pylon
Then, the design for the inverted Y-shaped pylon is presented. The design domain is shown in Fig. (6) and is discretized in a 30×80 mesh. The passive region is composed of the elements that cannot be optimized. It can either be void or solid from the beginning to the end. At the center of the top side, a loading of 5000 kN is applied. A volume fraction of 40% is required. As shown in Fig. (7), the optimized result is similar to the real pylon. The Rama VIII bridge in the Thailand is shown in Fig. (8). The resemblance between the optimized result and the Rama VIII bridge is noticeable. They both are typical Y-shaped pylons.
Original optimization model for inverted Y-shaped pylons.
Optimization result for inverted Y-shaped pylons.
The Rama VIII bridge.
If there are two loadings of 2500 kN each at the upper side, the design result will be different. The design domain is shown in Fig. (9), and the volume fraction is 35%. The optimized result is depicted in Fig. (10).
Design domain of suspension bridge tower.
Optimization result for suspension bridge tower.
For the sake of improving the accuracy of the conventional level set method and to assist engineers in design jobs with the level set method, a revised level set method for topology optimization is proposed using the trapezoidal method and topological sensitivity. Finite element analysis is utilized to discretize the structure. Based on result of the finite element analysis, the shape sensitivity and topological sensitivity are calculated to push and pull the boundary of material domain. In the conventional level set method, an initial guess of the number of internal holes is critical for its strong influence on the final optimized result. In the proposed method, the internal holes are nucleated spontaneously without an initial guess of the number of internal holes, and it manifests itself in the design of the cantilever, pylon and bridge tower. A flow chart of the algorithm implementing the proposed method is also provided. At the end, the design of the cantilever is presented to validate the proposed method. To demonstrate the utility, the proposed method is applied to the design of the inverted Y pylon and bridge tower. The optimized result for the inverted Y pylon design matches the practical designs very well.
The presented work only deals with the pylon and bridge tower design with simple loading case. The future work would extend the method to solve different optimization problems, such as dynamics issues. Also, the future work should address the issue of a structure that is meshed with irregular grids.
LIST OF ABBREVIATIONS
= a closed round sphere of radius ε
= an open round sphere of radius ε
= the boundary of
= the boundary of domain ∂D
= the boundary of which has a round sphere of radius ε
= the boundary of Ω
= the total computational domain for the optimization problem
= a positive number that approximates 0
= the strain tensor
= the external loading vector
= the non-linear constraint
= the Heaviside function
= an approximate Heaviside function
= the linear constraint
= the bandwidth of the approximate Heaviside
= the objective function
= the elemental stiffness matrix
= the global stiffness matrix
= the Lamé's first parameter
= the Lamé's second parameter
= the Lagrangian multipliers
= the velocity term defined by topological derivative
= the material domain
= a shape functional
= the level-set function for the domain
= the open bounded domain with a round sphere of radius
= the element densities field
= time, used in the evolution equation for the level-set function
= the stress tensor
= the topological derivative of
= the topological derivative of objective function
= the topological derivative of volume constraint
= the topological derivative of Lagrangian
= the displacement vector
= the elemental displacement vector
= the normal component of velocity
= the design velocity field
= the velocity term defined by shape derivative
= the coordinates of a point in the domain D
= the coordinates of the center of
CONFLICT OF INTEREST
The authors confirm that this article content has no conflict of interest.
This work was supported by the National Natural Science Foundation of China under Grant No. 51175197.
J.N. Richardson, S. Adriaenssens, P. Bouillard, and R.F. Coelho, "Multiobjective topology optimization of truss structures with kinematic stability repair", Struct. Multidiscipl. Optim., vol. 46, no. 4, pp. 513-532, 2012. [http://dx.doi.org/10.1007/s00158-012-0777-5]
K. Maute, and G.W. Reich, "Integrated multidisciplinary topology optimization approach to adaptive wing design", J. Aircr., vol. 43, no. 1, pp. 253-263, 2006. [http://dx.doi.org/10.2514/1.12802]
A.R. Díaz, and M.P. Bendsøe, "Shape optimization of structures for multiple loading conditions using a homogenization method", Struct. Optim., vol. 4, no. 1, pp. 17-22, 1992. [http://dx.doi.org/10.1007/BF01894077]
Z. Luo, N. Zhang, J.C. Ji, and T. Wu, "A meshfree level-set method for topological shape optimization of compliant multiphysics actuators", Comput. Methods Appl. Math., vol. 223–224, pp. 133-152, 2012.
M.Y. Wang, X. Wang, and D. Guo, "A level set method for structural topology optimization", Comput. Methods Appl. Math., vol. 192, no. 1, pp. 227-246, 2003.
M. Burger, "A framework for the construction of level set methods for shape optimization and reconstruction", Interface Free Bound, vol. 5, no. 3, pp. 301-330, 2003. [http://dx.doi.org/10.4171/IFB/81]
N.P. Van Dijk, G.H. Yoon, F. Van Keulen, and M. Langelaar, "A level-set based topology optimization using the element connectivity parameterization method", Struct. Multidiscipl. Optim., vol. 42, no. 2, pp. 269-282, 2010. [http://dx.doi.org/10.1007/s00158-010-0485-y]
N.P. Van Dijk, M. Langelaar, and F. Van Keulen, "Explicit level-set-based topology optimization using an exact heaviside function and consistent sensitivity analysis", Int. J. Numer. Methods Eng., vol. 91, no. 4, pp. 67-79, 2012. [http://dx.doi.org/10.1002/nme.4258]
F. Van Keulen, R.T. Haftka, and N.H. Kim, "Review of options for structural design sensitivity analysis. Part 1: Linear systems", Comput. Methods Appl. Math., vol. 194, no. 30–33, pp. 3213-3243, 2005.
V.J. Challis, "Multi-Property Topology Optimisation with the Level-Set Method,", Ph.D. thesis, The University of Queensland, Brisbane St Lucia, QLD, Australia, 2009.
G. Allaire, F. Jouve, and A.M. Toader, "Structural optimization using sensitivity analysis and a level-set method", J. Comput. Phys., vol. 194, no. 1, pp. 363-393, 2004. [http://dx.doi.org/10.1016/j.jcp.2003.09.032]