Three-Dimensional Dynamic Topology Optimization with Frequency Constraints Using Composite Exponential Function and ICM Method

The dynamic topology optimization of three-dimensional continuum structures subject to frequency constraints is investigated using Independent Continuous Mapping (ICM) design variable fields. The composite exponential function (CEF) is selected to be a filter function which recognizes the design variables and to implement the changing process of design variables from “discrete” to “continuous” and back to “discrete.” Explicit formulations of frequency constraints are given based on filter functions, firstorder Taylor series expansion. And an improved optimal model is formulated using CEF and the explicit frequency constraints. Dual sequential quadratic programming (DSQP) algorithm is used to solve the optimal model. The program is developed on the platform of MSC Patran & Nastran. Finally, numerical examples are given to demonstrate the validity and applicability of the proposed method.


Introduction
Structural topology optimization is to find optimal material layout within a given design space, for a given set of loads and boundary conditions such that the resulting layout meets a prescribed set of performance targets.The essence of topology optimization lies in searching for the optimum path of transferring loads; therefore the computational results of topology optimization are usually more attractive and more challenging than the results of cross-sectional and shape optimization.Although topology optimization is only in conceptual design phase in engineering, the design results significantly impact the performance of the final structure.
In the last decades, since the landmark paper of Bendsøe and Kikuchi [1], numerical methods for topology optimization of continuum structures have been developed quickly in application [2][3][4].Homogenization methods and SIMP (solid isotropic material with penalization) method are especially popular in continuum topology optimization, though their design variables are "fictions." In homogenization method [5,6], the design variables are microscale void and the optimal problem is defined by seeking the optimal porosity of porous medium using the optimality criteria.But sometimes the designs result in infinitesimal pores in the materials and the structure is not easy to manufacture.SIMP method [7][8][9] is employed to describe the relationship between the relative density and the material elastic modules of grey elements by a nonlinear function with continuous variables.These grey elements have intermediate densities by introducing a penalty coefficient.ESO method is another engineering approach which is investigated and used widely in recent years.The optimal topology is generated by deleting the group of elements with low strain energy from entire domain systematically [10][11][12].Later, a new development in ESO is the recent introduction approach named bidirectional evolutionary structural optimization (BESO) where elements are allowed to be added as well as removed [13,14].The validity and new-look of ESO type methods are stated by Huang and Xie [15][16][17], and a critical evaluation as well as comparison of SIMP and ESO method is discussed in detail by Rozvany [18].Other promising new methods include level Mathematical Problems in Engineering set method [19][20][21][22][23] and phase field method [24][25][26] and they are used widely in multiple material topology optimization.
Compared with static topology optimization, the study on dynamic topology optimization is more complicated and there are few investigations.Frequency topology optimization is of great importance in dynamic topology optimization and engineering fields.Topology optimization with respect to frequencies of structural vibration was first considered by Díaaz and Kikuchi [27], who studied the topology optimization of eigenvalues by using the homogenization method where reinforcement of a structure was optimized to maximize eigenvalues.Subsequently, many researches focus on expanding topology optimization in dynamic problems.Ma et al. [28,29] studied multi-eigenvalue optimization problems and frequency response optimization problems for vibrating plate structure.An optimization algorithm is derived to maximize a set of eigenvalues.Kosaka and Swan [30] investigated the optimization of the first five eigenfrequencies of plate structure and obtained the elastic properties by using a pure Reuss formulation.Krog and Olhoff [31] and Jensen and Pedersen [32] utilized a variable bound formulation that facilitates proper treatment of multiple frequencies.Pedersen [33] dealt with maximum fundamental frequency design of plates and included a technique to avoid spurious localized modes.Jensen and Pedersen [32] presented a method to maximize the separation of two adjacent frequencies in structures with two material components.Zhu and Zhang [34] emphasized on layout design which was related to optimization of boundary conditions and it was studied to maximize natural frequency of structures.In 2007, Du and Olhoff [35] introduced SIMP method for maximization of first eigenvalue and frequency gaps.Recently, Niu et al. [36] proposed a twoscale optimization method and found the optimal figurations of macrostructure-microstructure of cellular material with maximum structural fundamental frequency.Huang et al. [17] investigated the maximization of fundamental frequency of beam, plane, and three-dimensional block by applying a new bidirectional evolutionary structural optimization (BESO) method and dealt with localized modes by modifying the traditional penalization function of SIMP method.Xia et al. [37,38] presented a level set based shape and topology optimization method for maximizing the simple or repeated first eigenvalues of structure vibration.Further development on frequency topology optimization is seen in [3,[38][39][40][41].
However, in many researches concerning the model of dynamic topology optimization for continuum structure, frequency has been used as objective instead of constraints.In this paper, we extend our previous researches [42,43] primarily on static topology optimization issues of continuum structures to dynamic topology optimization field.Using structural weight as objective function, we build topology optimization model for continuum structure with frequency constraints.This method is beneficial to maintain the consistency of objective and constraint in cross-sectional optimization, shape optimization, and topology optimization.
Among the methods of mathematic optimization model solving, mathematical programming (MP) method is popular.Sequential quadratic programming (SQP) [44,45] in the MP method is widely used.SQP algorithm is a local quadratic approximation.A Lagrange function is needed to build for the global nonlinear optimization problem.However, a large amount of calculation is required if the number of optimization variables increases in the problem because of the repeated building of the Hessian of the Lagrange function.In this paper, dual sequential quadratic programming (DSQP) algorithm is employed to solve the topology optimization model constrained by frequency.In the process of solving optimal model, the dual theory is used to convert the constrained optimization model to one with reduced number of design variables, and the solving efficiency is greatly improved.
This paper is organized as follows.In Section 2, Independent Continuous Mapping (ICM) method is introduced.In Section 3, an improved frequency topology optimization model based on ICM method is built.Optimal algorithm to solve the mathematical optimization problem is given in Section 4. Two numerical simulations are presented in Section 5.In Section 6, conclusions are given.

Independent Continuous Mapping (ICM) Method
Independent Continuous Mapping (ICM) method is proposed by Sui [42] for skeleton and continuum structural topology optimization.It generalizes topological variables abstractly independent of the design variables such as sectional sizes, geometrical shape, density, or Young module of material.Using filter functions it maps the topological variables essentially with discrete values of zero and one to continuous independent topological variables with values belonging to [0, 1].
Based on the idea of ICM, the smooth model for structural topology optimization is established and it can be solved by the traditional algorithms in mathematical programming.Afterwards these continuous variables are discretized to 1 or 0 according to the criterion which the corresponding element is conserved or removed.The essential idea of the process is selecting the map from "discrete" to "continuous" and the inverse map from "continuous" to "discrete." For structural cross-section and shape optimization, natural frequency of structure is often taken as constraint.We denote   as the frequency of th order, and   ,   are the low and up bound of frequency, respectively.They satisfy the following inequality: (ii)   ≤   and  +1 ≥  +1 in nonfrequency band constraints.
Here  1 and  1 are natural and low bound of the first order frequency, respectively.For elastic structure, the usual relation between frequency f and eigenvalue  is  = (2) 2 .Therefore, the frequency constraints can apparently be transformed into eigenvalue constraints using the formula.Here we uniformly use () ≤  to generalize (i) and (ii) based on  = (2) 2 .
Thus, the model of continuum topology optimization with frequency constraints can be formulated as follows: where t, , and   denote the topological design variable vector, the total weight of structure, and the th element weight. and  are the th element and the th order frequency, respectively, and  and  represent the total number of constraints and elements.

Improved Model Based on ICM Method
3.1.Properties of Filter Function.In order to develop the model ICM method, we firstly investigate the essential part of ICM-the filter function.Its definition and choosing determine the establishment of optimization model and its solving, and further it will make great impact on the final performance of topology optimization.In order to map the topological variables from "discrete" to "continuous, " Sui [42] studied the filter function (  ) and summarized the following properties: (1) Continuity: filter function (  ) is continuous about topology variable   in the range of [0, 1]; (2) Differentiability: if   ∈ [0, 1], (  ) is continuously differential; (3) Monotonicity: (  ) is strictly monotonous; (4) Convexity: (  ) is strictly convex in [0, 1]; (5) Approximation: the approximation is given as follows: here  is an arbitrary small positive.However, the approximation property is not easy to control; it is a limit process.To achieve the approximation, we introduce the following formula: Consider ∫ 1 0 (  ) < ; here  is a given arbitrary small positive and 0 <  < 1/2.

Improved Filter Functions.
Several types of filter function are suggested in ICM method, among which Power Function (PF) is used frequently in application [46] and is as follows: Here   denotes th design variable. is a positive constant.
We introduce a new filter function-composite exponential function (CEF) to take the place of the old one and it is as follows: is a given positive constant and it is determined by numerical experiments with different problems.In Section 5, we compare the performance of the two types of filter function.

Improved Model of Frequency Topology Optimization.
Denote   (  ),   (  ), and   (  ) as filter functions for frequency topology optimization and they are given as follows: Here  0  , k 0  , and m 0  are the element weight, element stiffness matrix, and element mass matrix of original structure before the process of topology optimization, respectively.  , k  , and m  are the ones in the process of topology optimization, respectively.
Based on the above analysis, we build the model of topology optimization with frequency constraints as follows: In what follows, we face the problem of approximating the function of frequency constraints explicitly.
In the finite element analysis the dynamic behavior of a continuum structure can be represented by the following general eigenvalue problem: where K is the global stiffness matrix and M is the global mass matrix.  is the th eigenvalue and u  is the eigenvector corresponding to   .The eigenvalue   and the corresponding eigenvector u  are related to each other by Rayleigh quotient Since eigenvalue   is implicitly related with topology variable t, we use first-order Taylor series expansion for eigenvalue to express their relationship explicitly.
Take the reciprocal of stiffness filter function as design variables as follows: We have Therefore, the stiffness matrix filter function, mass matrix filter function, and weight filter function are given as follows: The global stiffness matrix K and the mass matrix M can be calculated by In view of ( 9) we have the derivative of   to design variable as follows: Substituting ( 13) to ( 14), we have where In (15),   = (1/2)u   k  u  and   = (1/2)  u   m  u  represent the strain energy and the kinetic energy of th element corresponding to the jth eigenmode, respectively.At this moment, the derivatives of eigenvalue with respect to all design variables can be obtained by subtracting the strain energy and kinetic energy for element mode from the results of modal analyses.
Using the first-order Taylor series expansion, the approximate expression of eigenvalue   (t) can be obtained: where superscript ] is the number at the ]th iteration.
If we define  as and further define then frequency constraints can be simplified by the following inequality: Thus ends the process of explicit approximation of the frequency constraints.

Solution of the Improved Topology Optimization Model
However, model ( 7) is a programming with nonlinear objective and linear constraints following the explicit process of frequency constraints.We use second-order Taylor series expansion to approximate the objective function and ignore the constant item.The model is transformed into the following quadratic programming model: As objective function varies with different filter functions, investigation of the different cases following different types of filter functions is necessary.Here we focus on PF and CEF.When PF is applied as the filter function, it is given as follows: where   ,   , and   are all arbitrarily integers greater than zero.In view of ( 22), we have   = 1/ Therefore the objective function ( 7) can be rewritten as where   = ( + 1) 0  /2(  ) +2 and   = −( + 2) 0  /(  ) +1 are undetermined parameters.
When CEF is applied as the filter function, it is given as follows: We have   = 1/  (  ) = ( 1/  − 1)/(   /  − 1), and therefore Similarly, the objective function in ( 7) can be expressed as where They are two undetermined parameters.The dual model of ( 21) is where z is the design variables in dual model: Here (x, z) is a Lagrange augmented function in (28).

Start
Construct the finite element model and input the optimal parameters Finite element analysis with MSC.Nastran.
Formulate the mathematical optimal model ( 7) Solve (21) with DSQP and update variables

No
Yes Then DSQP is employed to solve (28), and the precision of convergence is controlled to be 0.001.With the above analysis and solving of ( 28), the optimal topology structure with continuous design variables is obtained.

Numerical Simulation
Based on the above analysis, we outline the following computational procedure with a flow chart (Figure 1); it is used for maximizing the weight of structure constrained by frequency.
Following the above flow chart, we develop a program based on the platform of MSC Patran & Nastran to solve the improved topology optimal model.In order to demonstrate the effectiveness of the proposed method, two structures are presented as inputs.The first one is a cantilever with fundamental frequency constraint; the second one is a 3D beam with two concentrated masses by multiple frequency constraints.For the computation, the initial values of topology variables are all given as unit ( = 1); the lowest bounds of topology variables and the convergence precision values are 0.01 and 0.001, respectively.The solving topology configuration of the cantilever with different filter functions is given in Figure 3. Figures 3(a  half structure corresponding to (a)-(b).The performance of topology optimization with different filter functions is given in Table 1.The iterative curve of computation with different filter functions is described in Figure 4.  in the optimization process is presented in Figure 8 with different filter functions.

Conclusion
In this paper, an improved frequency topology optimization model of three-dimensional continuum structure is developed based on ICM method.CEF, a new filter function, is selected to recognize the design variables, as well as to implement much better the changing process of design variables from "discrete" to "continuous" and back to "discrete." Explicit formulations of frequency constraints are given based on filter functions and first-order Taylor series expansion and by extracting structural strain and structural kinetic energy from the results of structural modal analysis.An improved optimal model is formulated using CEF and the explicit frequency constraints.The program based on DSQP for solving the optimal model is developed on the platform of MSC Patran & Nastran.Finally, two examples of three-dimensional continuum structure show that clear and stable configurations can be obtained using ICM method.We find that configurations computed with DSQP combined PF and DSQP combined CEF are similar between one and the other in the case of fundamental frequency constraint.But we can find that DSQP combined CEF has the best performance for the optimization example from the point of view of optimal objective and iterative numbers.We also compute the first three order modal shapes of optimal structure with two different filter functions in multiple frequency constraints.The vibration modes from two filter functions are the same and the iteration curves of the first three order frequencies show that structural frequencies satisfy the frequency constraints.Therefore, all the computing results indicate that the proposed method is also feasible for multiple frequency constraints.

Figure 1 :
Figure 1: Flow chart of the iterative solution procedure.

5. 1 .
Example 1. Example 1 is a cantilever with size 100 × 50 × 50 mm 3 .It is fixed at one end as shown in Figure 2, and a concentrated mass Mc = 50 g is attached at the center of the other end.Young's modulus  = 100 GPa, Poisson's ratio  = 0.3, and mass density  = 1000 kg/m 3 .The structure is divided into 30 × 16 × 16 eight-node hexahedron elements.The 5010 Hz fundamental frequency of the cantilever is calculated by finite element analysis, and we set the frequency constraint for the design problem as  1 ≥ 3500 Hz.

Table 1 :
Performance of topology optimization with different algorithms and filter functions.

Table 2 :
Performance of topology optimization with different algorithms and filter functions.