Abstract

This paper shows how the reanalysis method can be utilized for speeding up the optimized process in topology optimization related to vibrating structures for simple and multiple eigenfrequencies. The block combined approximation with shifting (BCAS) method is used for reducing the computational effort included in repeated solution of eigenvalue problem, which will dominate a lot of the CPU time, especially for large problems. By utilizing Level 3 Basic Linear Algebra Subprograms (BLAS), the computation efficiency of the BCAS method is improved. For achieving an accurate optimal result, two indicators are presented to control the approximate reanalysis procedure. The effectiveness of the proposed method is demonstrated by three numerical examples.

1. Introduction

In classical topological optimization, the main purpose is concerned with the material distribution that represents the optimal layout in the admissible design domain [1]. The frequency optimization for vibrating structures is avoidance of resonance for the structure subjected to external excitation frequencies within a given interval [24]. For the optimal design procedure, two main subjects are needed to be treated carefully. The first one is the multiple eigenfrequencies problem, where the corresponding sensitivity analysis has to be calculated based on a mathematical perturbation analysis. The other is the localized modes in low-density areas.

Due to various changes in design, the frequency optimization for structures requires repeated solution of eigenvalue problem which involves repetitive and tremendous calculations. The aim of reanalysis methods is to efficiently analyze the structural responses for changes in design without full computation in optimization. Recently, some research efforts have been carried out for speeding up computational topology optimization procedures. For solving a large number of linear systems, an iterative minimum residual method for three-dimensional (3D) minimum compliance problem was proposed by Wang et al. [5], where parts of the search space corresponding to a selected design cycle were reused in the solution of the next design cycle. Based on recycled preconditioning, Amir [6] presented multigrid preconditioners generated at only certain design steps and reused in subsequent steps. Some other related approaches are based on utilizing Combined Approximations (CA) method. The main feature of CA method is the combination of efficient local approximations and accurate global approximations [7]. The CA method was originally developed for linear reanalysis. Later, it had been extended for eigenvalue (vibration) problem [8]. However, for CA method, the lower mode shapes can be accurately obtained, and only one eigenpair is calculated each time [9]. To address these shortcomings, some significant progress on vibration reanalysis has been developed, namely, modified CA method [10] and block combined approximations with shifting (BCAS) [11]. Based on optimality criteria schemes and mathematical programming, Leu and Huang [12] presented a reduced basis method for linear reanalysis. Amir et al. [13] integrated CA method into the framework of topology optimization, and the reanalysis method was used for solving the linear system of equilibrium equations. Bogomolny [14] put forward an efficient approach for topology optimization for free vibrations, and the procedure developed is based on the CA method which was utilized for repeated eigenvalue problem analysis. To accelerate the optimal process, Xu et al. [15] presented reanalysis-based GA (Genetic algorithm) optimization and its utilizations in the design of trusses. For speeding up the design process, Zuo et al. [16] presented the utilization of CA method in GA-based structural optimization for eigenvalue problems. Amir et al. [17] proposed an efficient robust topology procedure based on reanalysis techniques, and the CA method was used as a particular case of the preconditioned conjugates gradient procedure to solve the linear system. For topology optimization of structures under harmonic excitations, Zheng et al. [18] used the BCAS method for repeated eigenvalue problem analysis involved in the optimization procedure; no multiple eigenvalue phenomena is considered and the implementation of criteria for stopping the approximate reanalysis procedure and the objective function are different with this study.

In this paper, the BCAS method is utilized for topology optimization problems of vibrating structures. Compared with the our published work [18], the presented study focuses on maximum values of simple and multiple eigenfrequencies of vibrating structures, and the multiple eigenvalue phenomenon has been discussed in the optimization procedure. For increasing execution efficiency, Basic Linear Algebra Subprograms (BLAS) are applied. The BLAS include Level 1(perform scalar, vector and vector-vector operations), Level 2(perform matrix-vector operations), and Level 3 (perform matrix-matrix operations). It is worth noting that the BCAS method is based on matrix-matrix operations with Level 3 BLAS and hence the execution efficiency is improved [19].

The layout of the article is as follows. In Section 2, the formulations of optimization problem for eigenfrequency are reviewed, including the sensitivity analysis and the polynomial interpolation model. In Section 3, a brief description of the BCAS method for completeness is provided. In Section 4, some computational considerations are presented, and then three numerical examples are demonstrated for the accuracy and efficiency in Section 5. Finally, conclusions are presented.

2. Topology Optimization Formulation and Sensitivity Analysis

2.1. Topology Optimization Problems

The topology optimization problem for maximization of the eigenfrequency of vibrating structures can be stated as [2] subject towhere is the lower bound for the eigenfrequencies. and are the stiffness matrix and mass matrix, respectively. The eigenfrequency and the corresponding eigenvector can be obtained by solving (3). represents the number of the candidate eigenpairs, and it is suggested to be larger than that prevents the eigenfrequency exchange its order with the neighbouring eigenfrequency during the optimization process [2]. Let these eigenfrequencies be ordered so that and the corresponding eigenvectors satisfy constraint equation (4), where is Kronecker’s delta. In (5) and (6), the design variable is an element density and is used to prevent the mass and stiffness matrices from being singular. and denote the number of finite elements and the volume of element , respectively. is the prescribed structural volume.

With a view to avoiding the localized modes in the low-density areas, the Polynomial Interpolation Scheme (PIS) [20] is implemented in this paperwhere and are the stiffness matrix and mass matrix of element . and represent element stiffness matrix and mass matrix of solid element .

2.2. Sensitivities Analysis of Simple and Multiple Eigenfrequencies

In the following, the sensitivities of simple and multiple eigenfrequencies of vibrating structures with respect to the design variable can be derived. Note that multiple eigenfrequencies phenomenon may occur in some ways in structural optimization problems. One possibility is that an eigenfrequency is multiple at the beginning of the optimization procedure because of the structural symmetries, and the other possibility is that a simple eigenfrequency may also be multiple during the optimization procedure [2].

Firstly, consider the eigenfrequency is simple; i.e., , and then the corresponding eigenvector is unique. By directly differentiating the constraint equation (3) with respect to , we havePremultiplying (8) by and combining with (3) and (4), the sensitivity of a simple eigenvalue can be obtained bywhere the derivatives of the matrices and can be obtained from (7).

When the generalized eigenvalue problem in (3) has an -fold multiple eigenvalues , with and , that is to say, , corresponding eigenvectors are denoted by . In the case of , the calculation for the derivatives of eigenvalues is not straight forward. Based on a mathematical perturbation technique, a change of one arbitrarily selected design parameter is firstly considered, and corresponding stiffness and mass matrices will be changed; i.e., Then the multiple eigenvalues and corresponding eigenvectors for the perturbed design can be written aswhere and are unknown sensitivities of eigenvalue and corresponding eigenvector, respectively. Substituting (11) into (3), the first approximation can be obtained:Note that any linear combination of eigenvectors will satisfy (3), where are unknown coefficients. Premultiplying (12) by hasand we haveFinally, the sensitivities of multiple eigenvalues are determined byFor more detailed explanation for sensitivity analysis of the multiple eigenvalues, readers are referred to [21].

To avoid the checkerboards pattern and mesh-dependencies phenomena, the mesh-independency filter technique [22] is adopted to modified the sensitivity asand the weight factor in (16) is defined aswhere the operator is defined as the distance between the center of element and the center of element , and is the prescribed filter radius.

3. Vibration Reanalysis by BCAS

In this section, the approximate reanalysis is implemented following the block combined approximation method with shifting (BCAS) [11]. Here, the formulation of freely vibration reanalysis by BCAS is briefly described for completeness. Consider the generally eigenvalue problem resulting from a finite element (FE) discretization of a “reference” design or the initial design corresponding to a certain iterative circle of the optimization procedurewhere and are the stiffness and mass matrices, respectively. The lowest eigenvalues and corresponding can be calculated via subspace iteration with shifting [23] or shifted Lanczos algorithm [24]. In addition, a series of factorizations of the shifted stiffness matrices is involved in the solution procedure; is the number of these triangular factorizations.

According to the BCAS, a modified eigenvalue problem is introduced as follows:where , , and and represent the changes of stiffness and mass matrices due to the changes in design, and the corresponding generalized eigenvalue problem with shifting can be given as where . The two eigenvalue problems in (19) and (20) have the same eigenvectors. It can be derived from (20) thatand (21) can be expanded aswhereSubstituting for and dropping in (22), the basis vectors in the BCAS method are evaluated by

As the decomposed form of is usually given from the initial design, the calculation of only contains forward and backward substitutions. The stiffness matrix and mass matrix are condensed asThe modified analysis equation (19) can be approximated by the following reduced eigenvalue problemwhere and are the matrices of eigenvalues and the corresponding eigenvectors. The matrix of eigenvalues and corresponding matrix of eigenvectors for the modified analysis equation (19) can be evaluated byFor BCAS method, the approximate solutions can be high accuracy achieved with . By using the information of the eigenvalues for the initial design, the shifted frequency in (20) can be chosen asWhen the number of required eigenpairs is small, we usually set . Note that the reanalysis method is based on matrix-matrix operations and can provide higher performance by using the Level-3 BLAS. More detailed explanation can be found in [11].

4. Computational Considerations

The major aim of utilizing approximate reanalysis is to obtain an accurate result efficiently. However, large changes in stiffness matrix may result in the large error. Hence the key for efficiently achieving an accurate result is choosing the right time to stop a series of approximate reanalysis and implement an exact solution [13]. In this paper, two indicators are adopted to control the approximate optimization procedure: the absolute error between and is defined asand the Modal Assurance Criterion (MAC) value is defined by where , are the design variable vectors at step and the one corresponding to the previous exact solution, respectively.

The convergence criteria [25] of the optimization procedure are set to where is the objective value of step. Meanwhile, is used to illustrate the accuracy of the result, where and denote the final design variable vectors calculated by the proposed procedure and exact method, respectively. The implementation of the proposed optimization procedure is shown in Figure 1.

5. Numerical Implementation

In this section, three numerical examples are given to demonstrate the efficiency of the proposed method. The stiffness and mass matrices are stored in the Harwell-Boeing (HB) compressed format. The subspace iteration (SI) method with shifting is used to solve the eigenvalue problem in exact method. In these numerical examples, must be chosen to be sufficiently large to include all members of a possible -fold eigenfrequency, such as or . The Globally Convergent Method of Moving Asymptotes (GCMMA) algorithm [26] is used for the solution of the optimization problem. The related computations are based on Intel Visual Fortran XE 2015 with Math Kernel Library 11.2, and efficient Direct Sparse Solver (DSS) Interface Routines-sparse storage, factorization, and sparse solver are used. The computations of all examples are completed on the platform: Intel CPU Xeon E5-2623 v4, 32G RAM.

5.1. Maximization of the First Eigenfrequency of a 2D Cantilever Beam

For verifying the efficiency of the proposed method, we use the numerical example presented in [27] as the first example. Figure 2 shows a rectangular domain of , and its thickness is . The structure is clamped at the left side. The design domain is divided into four-node quadratic elements and total DOFs (degrees of freedom) are about 6560. Young’s modulus is assumed , Poisson’s ratio , and the mass density . A concentrated nonstructural masse is applied in the corners of the right edge.

The objective function is to maximize the first eigenfrequency and the volume of available material is constrained to be less than 40%. The converged results via the proposed method and exact method are given in Figure 3. It indicates that there is no significant difference of the optimized configurations of the structure based on the two methods. The convergence curves of the objective function are given in Figure 4. The more details of optimization are listed in Table 1. The exact method reaches objective value of 77.09 rad/s while the proposed method reaches 76.35 rad/s (1% error), and the MAC value for the design variables based on exact method and the proposed method is 0.982. The proposed method can thus achieve very excellent approximation to the exact result. In addition, CPU time for the topology optimization based on the BCAS method could be decreased to 35.3% as the number of the SI method reduced. Note that the final design is similar to the one presented in [18], but the objective function of this example is different.

5.2. Maximization of the First Eigenfrequency of a 2D Rectangle Clamped Beam

The beam of dimension is clamped on both ends as shown in Figure 5, and its thickness is . It is divided into four-node quadratic elements with 22878 DOFs. Young’s module of material is set to , Poisson’s ratio , mass density . A concentrated nonstructural mass is placed at the center.

The objective function is to maximize the first eigenfrequency of the structure with 40% material of the design domain. Similar final optimization results achieved by the proposed method and exact method are shown in Figure 6. Figure 7 shows the convergence curves of the objective function based on the proposed method and the exact method. Performance results of topology optimization with exact method and the proposed method are shown in Table 2. From Table 2, it can be observed that the objective value obtained by exact method is 126.91 rad/s, and the one obtained by the proposed method is 123.56 rad/s (2.6% error). Meanwhile, the MAC value between the final design variables by the proposed method and exact method is 0.975. Moreover, it can be found that the topology optimization based on BCAS can reduce 34% CPU time.

5.3. Maximization of the Second Eigenfrequency of 3D Clamped Square Plate

The structure has the domain of and is clamped at four edges as shown in Figure 8. It is divided into 3D 8-node solid elements with 13689 DOFs. Young’s module of material is , Poisson’s ratio is , and the mass density is . An additional artificial mass is applied at the center of the structure.

Here, the objective function is to maximize the second eigenfrequency of the structure, and the volume fraction is constrained to be less than 50%. The calculation of the sensitivity numbers for multiple eigenfrequencies should be taken into consideration because of structural symmetry. The similar final topology designs obtained by exact method and the proposed method are introduced in Figures 9(a) and 9(b). Figure 10(a) shows the first three eigenfrequencies as a function of design iterations. It can be observed that the multiple eigenfrequencies phenomenon appears at the beginning of iteration; i.e., the second eigenfrequency is superposed with the third eigenfrequency with . Figure 10(b) gives the comparison of iteration histories for the objective function based on the proposed method and exact method. Details of computational result are listed in Table 3. The final objective value obtained by the proposed method is 175.95 rad/s (1.0% error), and the MAC value is 0.996 which reflects a good correlation between the final design variables obtained by the proposed method and exact method. Moreover, it can be found that, compared with the CPU run-time of exact method, the one of the proposed method is decreased by 32%.

6. Conclusions

In this paper, the BCAS method for freely vibration reanalysis, which can be used to compute several or more eigenpairs of the modified structure, has been applied to topology optimization of vibrating structures for simple and multiple eigenfrequencies. Two indicators are adopted to control the approximate reanalysis procedure. The multiple eigenfrequencies phenomenon is considered during the optimization procedure. The BCAS method is based on matrix-matrix operations with Level 3 BLAS and hence the execution efficiency is improved. Numerical examples have shown successful implementation of the topology optimization based on BCAS method.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to thank Svanberg K. for providing the Matlab code of GCMMA optimizer. The work was supported by the National Natural Science Foundation of China (Grant no. 41602369), the Science and Technology Developing Plan Project of Jilin Province (Grant no. 20160520021JH), and Key Laboratory of Ministry of Land and Resources on Complicated Conditions Drilling Technology (Grant no. KLDET2017013).