Plate/shell topological optimization subjected to linear buckling constraints by adopting composite exponential filtering function

In this paper, a model of topology optimization with linear buckling constraints is established based on an independent and continuous mapping method to minimize the plate/shell structure weight. A composite exponential function (CEF) is selected as filtering functions for element weight, the element stiffness matrix and the element geometric stiffness matrix, which recognize the design variables, and to implement the changing process of design variables from “discrete” to “continuous” and back to “discrete”. The buckling constraints are approximated as explicit formulations based on the Taylor expansion and the filtering function. The optimization model is transformed to dual programming and solved by the dual sequence quadratic programming algorithm. Finally, three numerical examples with power function and CEF as filter function are analyzed and discussed to demonstrate the feasibility and efficiency 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 B Hong-Ling Ye yehongl@bjut.edu.cn 1 College of Mechanical Engineering and Applied Electronics Technology, Beijing University of Technology, Beijing 100124, China 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 crosssectional and shape optimization. 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]. The classical methods include the homogenization method [5,6], the variable density method (including solid isotropic material with penalization model (SIMP) and rational approximation of material properties (RAMP) interpolation model) [7][8][9][10], evolutionary structural optimization (ESO) [11][12][13], level set method [14][15][16], and so on.
The plate/shell structure is popular for lightweight constructions in national defense and civil industries. However, it is shown from both research literature and industrial applications that plate/shell structures are prone to buckle. As buckling affects the security of the whole structure, it is necessary to address the structural stability during structure design. Buckling topology optimization of plate/shell is to find optimal material layout of plate/shell structure that meets the buckling requirements. Although buckling topology optimization is only in the phase of conceptual design in engineering, the optimal results will impact the performance of the final structure significantly. Compared with static topology optimization, buckling topology optimization is more complicated, and there are few investigations up till now. In 1995, linear buckling topology optimization of two-dimensional structures had been studied by Neves et al. [17], the optimization results lay the foundations for the non-linear buckling optimization. Meanwhile, Seo [18] studied the topology optimization of inner-wall stiffener of cylindrical containers. The reciprocal of critical buckling load was adopted as an objective function, and the total mass of stiffener was constrained to a prescribed value. Later, Neves et al. [19] presented a two-scale asymptotic method for the linearized elastic stability analysis. Topology optimization of the periodic microstructures is carried out based on the local buckling instabilities in the periodic boundary conditions (PBC). Combining the linearized elastic buckling model with the inverse homogenization and an eigenvalue buckling analysis with Floquet-bloch wave theory [20], the minimum critical buckling strain is obtained and maximized with the PBC having a constant volume fraction of materials. In 2002, Ramm et al. [21,22] constructed the topology optimization model with linear buckling constraints based on the SIMP method to study the influence of geometrical nonlinear behavior to topology optimization design. In 2004, the shell's topological optimization under linear buckling response using SIMP power-law penalization of stiffness was given to achieve the discrete topology [23]. In 2009, Lund [24][25][26] studied the buckling topology optimization of laminated multi-material composite shell structures by introducing interpolation functions, which is from the SIMP approach. In 2012, Browne [27] studied the method of solving the large-scale quadratic programming problem, and the method is applied to the topology optimization problem using compliance and buckling as constraints with the minimum structural weight as objective. In 2013, Lindgaard [28] studied the static geometry nonlinear structure topology optimization of instability to maximize the buckling critical load.
Up till now, different optimization methods have been used to solve the buckling structure topology optimization problems, however, there is no uniform valid method to deal with the buckling topology optimization of plate/shell. This is so since building the buckling topology optimization model is more complex and difficult than static topology optimization, and calculations for sensitive analysis are enormous.
In this paper, we investigate buckling topology optimization based on independent, continuous and mapping (ICM) method, proposed by Sui [29] for skeleton and continuum structural topology optimization in 1996. The topological variables are independent of design variables such as sectional sizes, geometrical shape, density or Young's modulus of material. Filter functions are used to map the changing process of topological design variables from "discrete" to "continuous" and back to "discrete". The smooth model for structural topology optimization is established and solved by the traditional algorithms in mathematical programming. The ICM method has been mainly used to study static and dynamic topology optimization [30][31][32]. We extend this method and do in depth research for the buckling problem. A model of topology optimization for the lightest plate/shell structures with the critical buckling load constraints is constructed. Usually, a power function (PF) is selected as the filter function in the past, and we select a composite exponential function (CEF) as the filtering function to complete the changing process of design variables. The optimal results with two different filter functions are compared by numerical examples.
This paper is organized as follows. In Sect. 2, a buckling topology optimization model of plate/shell structure based on the ICM method is established. In Sect. 3, a strategy for solving the buckling optimization model is introduced. An optimal algorithm to solve the mathematical optimization problem is given. In Sect. 4, the program flow of the optimization algorithm is charted. Three numerical simulations are presented in Sect. 5. In Sect. 6, conclusions are given.

Linear buckling analysis of plate/shell structures
Structural buckling widely exists in practical engineering structure. Buckling is a mathematical instability, leading to a failure mode. As the applied load is increased on a structure by a small amount beyond the critical load, the structure deforms into a buckled configuration. Further load will cause significant and somewhat unpredictable deformations, possibly leading to complete loss of the structural load-carrying capacity. The interpretation of this result is that for P < P cr j , the structure remains stable. For P > P cr j , the structure is unstable and buckles. P cr j is the critical load for buckling. Usually, once the form of structure is established, its buckling will have a variety of modes and multiple critical loads. The structure will not work before the mode reaches higherorder buckling mode, so we just care about the first-order critical load of buckling mode. In this paper, the linear elastic and pre-buckling of continuum structure is considered. Assuming the structure to be perfect with no geometrical imperfections, stresses are proportional to the loads, i.e., stress stiffness depends linearly on the load, displacements at the buckling configuration are small, and the load is independent of the displacements. The linear buckling problem can be represented as [28,29] where K and G denote the structural stiffness matrix and geometric stiffness matrix, respectively, λ j is the j-th eigenvalue, i.e., buckling critical load factor and u j describes the corresponding eigenvector, j denotes the j-th order of the modal.

Description of the filter function based on the ICM method
Filter function is the key strategy of the ICM method. It identifies the corresponding element or subdomains of geometric quantity or physical quantity, such as the element weight, the element stiffness matrix, and the element geometric stiffness matrix. Discrete design variables can be mapped to continuous variables by filter function and inversed back to discrete variables. For the buckling topology optimization model, we define element weight, the element stiffness matrix, and the element geometric stiffness matrix as follows where t i is the topology variable value of the i-th element. w i , k i , g i denotes the element weight, stiffness matrix, and geometric stiffness matrix of i-th element in the optimal process, respectively. And w 0 i , k 0 i , g 0 i represent the initial element weight, stiffness matrix, and geometric stiffness matrix of the i-th element, respectively. f w (t i ), f k (t i ), and f g (t i ) are the filter functions of the element weight, element stiffness matrix, and the element geometric stiffness matrix.
In addition, the element weight, element stiffness matrix, element geometry stiffness matrix, and element quality matrix are changed by taking advantage of filter functions. These physical quantities of every element change a lot when the structural topology changes, and then the filter functions in the formulation can lead to convergence. Furthermore, filter functions influence the speed of convergence and the stability of the solution of the optimal process.
Several types of filter function are suggested in the ICM method [33]. Among which, the PF is used frequently as follows Here, α is a given positive constant. Now, we introduce a new filter function-CEF to take the place of PF, and it is as follows where γ is a given positive constant. In Sect. 5, the performances of the two types of filter function are compared. From Eqs.
(2), (3), and (4), the specific expressions of PF and CEF in the model of buckling topology optimization are given : It should be pointed out that these parameters of filter functions can be determined by using the least squares method or numerical experiments, see Refs. [29][30][31].

Mathematical model of buckling topology optimization
Based on the ICM method, the optimal model to minimize the structural weight subjected to the linear buckling constraints is as follows where t denotes the vector of topological design variables, W is the structural weight, and w i is the element weight of structure, P cr j presents the critical buckling load, P cr j is the lower limit buckling critical load, and J and N denote the total number of the buckling modes and number of elements.
The relationship between critical load and external load P can be expressed as Then the buckling critical load factor is used as constraints in the optimal model. The buckling topology optimal model (6) can be transformed as follows where λ j is the lower limit buckling critical load factor. In order to solve the optimal model, the reciprocal of filter function with stiffness matrix is used as a design variable as follows Therefore, the topological design variable is expressed as Then Eq. (2) can be transformed into Eq. (11) With the introduction of filtering functions and the reciprocal of filter function of stiffness matrix, the optimization model (8) is written as 3 Strategy for solving the buckling optimization model

Design sensitivity analysis
To estimate the design sensitivity, we have to consider the derivative of the eigenvalue in Eq. (1). The eigenvalues can be expressed using the Rayleigh quotient: The derivative of the eigenvalue is given as follows β i can be obtained according to different type of filter function. When PF is selected as filter function, And CEF acts as filter function, Therefore, Eq. (14) is deduced as Here, U i j = 0.5u T j k i u j and V i j = 0.5u T j g i u j represent the strain energy and geometric strain energy for the i-th element in j-th mode, respectively. V j = u T j Gu j is the structural geometric strain energy in the j-th mode.

Explicit approximation of buckling constraints
As the constraint is implicit about design variables, we make it explicit by using the first order Taylor expansion: Here, the superscript v is the number of iterations. Take Eq. (20) into Eq. (21), we have where

Then the buckling constraints in model (12) can be expressed as
We set So the buckling constraints can be simplified to the following form

The standardization of the objective function
In order to obtain an explicit objective, the second-order Taylor expansion is used. When PF is selected as a filter function, the structural weight can be written Here, When CEF acts as filter function, the structural weight can be written where, Therefore, the standard quadratic programming model for Eq. (12) can be obtained as follows

Solution of the optimization model
As the number of design variables is much bigger than that of the constraints, we deduce the dual model to decrease the number of design variables as follows in order to reduce the amount of calculation.
where z is the design variable vector of the dual model, Φ (z) is the objective function, and Φ(z) = min (L(x, z)), In this paper, the convergence criterion is chosen as follows W (v+1) and W (v) are the current iteration and the previous iteration of structural weight. ε is the convergence precision, ε = 0.001.

Discrete degree of topological design variables
In order to measure the discrete degree of topological design variables, we use M nd [34] as a criterion, and it is given where T i is the topological variable value for the i-th element and N is the total number of the elements. Following Eq. (29), if the topological variable value is 0 or 1, then M nd is 0; if the topological variable values is 0.5, then M nd is 1. The closer the topological variable value to 0 or 1, the smaller is the value of M nd and the better the optimal result.

Program flow of optimization algorithm
The numerical procedures are developed by the PCL toolkit in the MSC. Patran software platform. We use MSC.Nastran to analyze the numerical solution of Eq. (1). The corresponding program flow as shown below Step 1 Build finite element model by using MSC. Patran; Step 2 Input initial optimal parameters and set up optimal model; Step 3 Make buckling analysis by using MSC. Nastran; Step 4 Calculate and extract the critical buckling factor and strain energy; Step 5 Input parameters of the optimal algorithm; Step 6 Solve the dual optimization model (27) by the dual sequence quadratic programming (DSQP); Step 7 Judge convergence of the optimal results. If the structural weight satisfies the formula (28), then the program is terminated. Otherwise, update design variables x and topology design variables t, then go to step3.

Numerical examples
In this section, three examples of topology optimization of single material plate/shell structures are given. All the material is isotropic with Young's modulus E = 68890 MPa, Poisson's ratio μ = 0.3. In the initial design, the available material is uniformly distributed over the admissible design domain. The structures are meshed by four-node 2D plate/shell finite elements. The specific boundary condition and force form are shown in the following three examples. Fig. 1, the base structure is a plane elastic body with size 520 × 260 × 2 mm 3 , and mass density ρ = 1.0 × 10 −3 kg/mm 3 . The distributed compression load at the top edge is 100 N/mm. Two corners of the bottom edge are fixed. The buckling constraint value is 100 N in Ref. [34]. The region which including the above two layers of the unit is a non-design and should be maintained. The base structure's weight is 0.73 kg.

Example 1 As shown in
The topology configuration of the structure with CEF filter function is given in Fig. 2. It is similar to the topology configuration with PF as in Ref. [34] as shown in Fig. 3. The iterative history curve of buckling load and structural weight are shown in Figs. 4 and 5. The optimal structural weight  with CEF is 104.808 kg and the iterative number is 36, as the optimal structural weight with PF is 115.756 kg and the iterative number is 51. From the point view of structural weight and iteration, the optimal results with CEF is better that of PF. Fig. 6, the base structure is a plane elastic body and the mass density is ρ = 2.7 × 10 3 kg/m 3 . The forces P = 15000 N are located on the midpoint of the top and bottom boundary. Four corners of the structure are fixed. The base structure's weight is 0.73 kg.

Example 2 As shown in
After finite element analysis, the first-order buckling load factor of the structure is λ 1 = 0.0533. The critical buckling load is 1600 N, and the buckling load 1300 N is defined as constraint value.
The topology configurations of the structure with different filter functions before and after discretion are given in Fig. 7. In addition, the first-order buckling modal of optimal structure is computed as shown in Fig. 8. The iterative history of buckling load and structural weight with different filter functions are depicted in Figs. 9 and 10. We can see that the critical buckling loads of the optimal structures satisfy the buckling constraint. From Fig. 10, we can see a clear differ- The first-order buckling modal of optimal structure with different filter function. a PF. b CEF ence with PF and CEF in the structural weight and iterative number. The optimal structural weight with CEF is lighter and the number of iterations is less than that of PF.
The distributions of topological design variables are listed in Table 1. The discretization of topological design variables is evaluated by using M nd . We can find that M nd with PF and CEF are 11.60 % and 7.25 %. Therefore, the optimal result with CEF is better than that of PF from the view of discretization of topological design variables. Fig. 11, the base structure is a part of the cylindrical shell, the generatrix is 260 mm, arc length is 520 mm, and the radius is 5000 mm. The force P = 15000 N is located on the center of the cylindrical shell along the radial direction. After finite element analysis, the first-order   The intermediate results and optimal topology configuration of the structure with different filter functions are indicated in Fig. 12. The iterative history curve of buckling load and structural weight with PF and CEF is given in Figs. 13 and 14. The performances of topological optimization with different filter functions are given in Table 2. From Fig. 13, we get that the optimal structure with PF and CEF all satisfy the buckling constraint. The optimal structural weight with CEF is lighter than that of PF as shown in Fig. 14.

Example 3 As shown in
From the above three examples, we can see that the objective (weight) with CEF is apparent lighter than that of PF. We can also find that the optimal result with CEF has the best performance from the point of view of iterative number. The distribution of optimal topological values show that the percentages of M nd with CEF is lower than that of PF, so the CEF filter function has the best performance from the viewpoint of discreteness.

Conclusion
In this paper, a buckling topological optimal model of plate/shell structure is established based on the ICM method.  CEF is selected as a filter function to recognize the design variables, as well as to implement the changing process of design variables from "discrete" to "continuous" and back to "discrete". Explicit formulations of buckling constraints are given based on two different filter functions, first-order Taylor series expansion by extracting structural strain and structural kinetic energy from the results of structural modal analysis. The program based on DSQP for solving the optimal model is developed on the platform of MSC. Patran & Nastran.
Three numerical examples of continuum structure show that clear and stable configurations can be obtained by using the ICM method. We also find that configurations computed with PF and CEF are similar. But we can see that the optimal result with CEF has the better performance from the point of view of optimal objective, iterative numbers, and discrete degree.