An Efficient Reanalysis Method for Topological Optimization of Vibrating Continuum Structures for Simple and Multiple Eigenfrequencies

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 andmultiple 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 BCASmethod 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.


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 [2][3][4].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.

Topology Optimization Formulation and Sensitivity Analysis
. .Topology Optimization Problems.The topology optimization problem for maximization of the ℎ eigenfrequency of vibrating structures can be stated as [2] max  1 ,...,

{𝛽}
(1) subject to where  is the lower bound for the ℎ eigenfrequencies.K and M 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 0 <  1 ≤  2 ≤ ⋅ ⋅ ⋅ ≤   and the corresponding eigenvectors [ 1 ,  2 , . . .,   ] satisfy constraint equation ( 4), where   is Kronecker's delta.In ( 5) and ( 6), the design variable   ( = 1, . . .,   ) is an element density and  ≤ 0.001 is used to prevent the mass and stiffness matrices from being singular.  and V  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 lowdensity areas, the Polynomial Interpolation Scheme (PIS) [20] is implemented in this paper where k  and m  are the stiffness matrix and mass matrix of element .k 0 and m 0 represent element stiffness matrix and mass matrix of solid element .
. .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   ( = 1, . . .,   ) 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   (  =  2  ) is simple; i.e.,  −1 <   <  +1 , and then the corresponding eigenvector   is unique.By directly differentiating the constraint equation ( 3) with respect to   , we have Premultiplying ( 8) by    and combining with (3) and ( 4), the sensitivity of a simple eigenvalue   can be obtained by where the derivatives of the matrices K and M can be obtained from (7).
When the generalized eigenvalue problem in (3) has an -fold multiple eigenvalues   =  2  ,  = +1, +2, . . ., +, with  ≥ 0 and  +  ≤ , that is to say,  +1 =  +2 = ⋅ ⋅ ⋅ =  + = , corresponding eigenvectors are denoted by  +1 ,  +2 , . . .,  + .In the case of  > 1, 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 as where   and k  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) For more detailed explanation for sensitivity analysis of the multiple eigenvalues, readers are referred to [21].
To avoid the checkerboards pattern and meshdependencies phenomena, the mesh-independency filter technique [22] is adopted to modified the sensitivity   /  as and the weight factor   in ( 16) is defined as where the operator dist(, ) is defined as the distance between the center of element  and the center of element , and  min is the prescribed filter radius.

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 procedure where K and M are the  ×  stiffness and mass matrices, respectively.The  lowest eigenvalues Λ = diag( 2 1 ,  2 2 , . . .,  2  ) and corresponding Φ = [ 1 ,  2 , . . .,   ] can be calculated via subspace iteration with shifting [23] or shifted Lanczos algorithm [24].In addition, a series of LDL T factorizations of the shifted stiffness matrices K−  M( = 1, . . ., ,   ≥ 0) 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 K  = K + ΔK, M  = M + ΔM, and ΔK and ΔM 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 Λ  = Λ  −   I.The two eigenvalue problems in (19) and (20) have the same eigenvectors.It can be derived from (20) that Mathematical Problems in Engineering and ( 21) can be expanded as where Substituting Φ for Φ  and dropping Λ  in (22), the  ×  basis vectors R  = [R 1 , R 2 , . . ., R  ] in the BCAS method are evaluated by As the decomposed form of K −   M is usually given from the initial design, the calculation of R  only contains forward and backward substitutions.The stiffness matrix K  and mass matrix M  are condensed as The modified analysis equation ( 19) can be approximated by the following reduced eigenvalue problem where Λ  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 by For BCAS method, the approximate solutions can be high accuracy achieved with  = 2.By using the information of the eigenvalues for the initial design, the shifted frequency   in (20) When the number of required eigenpairs is small, we usually set   = 0. Note that the reanalysis method is based on matrixmatrix operations and can provide higher performance by using the Level-3 BLAS.More detailed explanation can be found in [11].

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 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.

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  =  + 2 or  =  + 3. The Globally Convergent Method of Moving Asymptotes (GCMMA) algorithm [26]  . .Maximization of the First Eigenfrequency of a D 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 8 × 4, and its thickness is 1.The structure is clamped at the left side.The design domain is divided into 80 × 40 four-node quadratic elements and total DOFs (degrees of freedom) are about 6560.Young's modulus is assumed  = 1.0 × 10 10 , Poisson's ratio ] = 0.3, and the mass density  = 1000g/ 3 .A concentrated nonstructural masse  = 16000g is applied in the corners of the right edge. 2. Detect possible N -fold multiple eigenvalues  2 j .

4.
Solve the generalized eigen-problem of the updated structure for the eigenfrequencies and modes by Finite Element analysis: if method; if not, perform the SI method.
No  e converged?Yes perform the postprocessor result Stop 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.

. . Maximization of the First Eigenfrequency of a D Rectangle
Clamped Beam.The beam of dimension 14 × 2 is clamped on both ends as shown in Figure 5, and its thickness is 0.01.It is divided into 280 × 40 four-node quadratic elements with 22878 DOFs.Young's module of material is set to  = 2.5 × 10 10 , Poisson's ratio ] = 0.3, mass density  = 2.5 × 10 3 g/ 3 .A concentrated nonstructural mass  = 125g 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   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 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%.

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.

Figure 1 :
Figure 1: The flow chart of the proposed optimization procedure.

Figure 3 :
Figure 3: The converged results for a 2D cantilever beam: (a) exact method; (b) proposed method.

Figure 4 :Figure 5 :
Figure 4: Convergence curve for the objective function of a 2D cantilever beam: proposed method versus exact method.

Figure 6 :
Figure 6: The converged results for a 2D rectangle clamped beam: (a) exact method; (b) proposed method.

Figure 7 :Figure 8 :
Figure 7: Convergence curve for the objective function of a 2D rectangle clamped beam: proposed method versus exact method.

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  2 =  3 = 65.37/.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
adopted to control the approximate optimization procedure: the absolute error between   and    is defined as        −    = [ ,1 , . . .,  ,  ],    = [  ,1 , . . .,   ,  ] are the design variable vectors at ℎ step and the one corresponding to the previous exact solution, respectively.
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.

Table 1 :
Results of maximization the first natural frequency of the 2D cantilever beam.

Table 2 :
Results of maximization the first natural frequency of the 2D rectangle clamped beam. = 1.0 × 10 11 , Poisson's ratio is ] = 0.3, and the mass density is  = 7800g/ 3 .An additional artificial mass  = 3.12 × 10 5 g is applied at the center of the structure.
. .Maximization of the Second Eigenfrequency of D Clamped Square Plate.The structure has the domain of 20 × 20 × 1 and is clamped at four edges as shown in Figure 8.It is divided into 40 × 40 × 2 3D 8-node solid elements with 13689 DOFs.Young's module of material is

Table 3 :
Results of maximization the second natural frequency of the 3D clamped square plate.