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

Authors: Youhong Sun; Xuqi Zhao; Yongping Yu; Shaopeng Zheng
Journal: Mathematical Problems in Engineering
 (2018)

---

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


---

## Body

## 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 [2–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 Section2, 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 theith eigenfrequency of vibrating structures can be stated as [2] (1)maxρ1,…,ρNEβsubject to(2)β-ωl2≤0,l=i,i+1,…,J(3)Kφj=ωj2Mφj,j=1,…,J,(4)φjTMφk=δjk,j,k=1,…,J,,(5)∑e=1NEρeve≤V-,,(6)0≤ρ_≤ρe≤1,e=1,…,NE..where β is the lower bound for the ith eigenfrequencies. K and M are n×n the stiffness matrix and mass matrix, respectively. The ith eigenfrequency ωi and the corresponding eigenvector φi can be obtained by solving (3). J represents the number of the candidate eigenpairs, and it is suggested to be larger than i that prevents the ith eigenfrequency exchange its order with the neighbouring eigenfrequency during the optimization process [2]. Let these eigenfrequencies be ordered so that 0<ω1≤ω2≤⋯≤ωJ and the corresponding eigenvectors [φ1,φ2,…,φJ] satisfy constraint equation (4), where δjk is Kronecker’s delta. In (5) and (6), the design variable ρe(e=1,…,NE) is an element density and ρ_≤0.001 is used to prevent the mass and stiffness matrices from being singular. NE and ve denote the number of finite elements and the volume of element e, respectively. V¯ 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 paper(7)ke=15ρe5+ρe16ke0me=ρeme0where ke and me are the stiffness matrix and mass matrix of element e. ke0 and me0 represent element stiffness matrix and mass matrix of solid element e.


### 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ρe(e=1,…,NE) 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 thejth eigenfrequency ωj(λj=ωj2) is simple; i.e., ωj-1<ωj<ωj+1, and then the corresponding eigenvector φj is unique. By directly differentiating the constraint equation (3) with respect to ρe, we have(8)K-λjM∂φj∂ρe+∂K∂ρe-λj∂M∂ρe-∂λj∂ρeMφj=0,e=1,…NEPremultiplying (8) by φjT and combining with (3) and (4), the sensitivity of a simple eigenvalue λj can be obtained by(9)∂λj∂ρe=φjT∂K∂ρe-λj∂M∂ρeφj,e=1,…,NEwhere the derivatives of the matrices K and M can be obtained from (7).When the generalized eigenvalue problem in (3) has an N--fold multiple eigenvalues λj=ωj2,j=s+1,s+2,…,s+N-, with s≥0 and s+N-≤J, that is to say, λs+1=λs+2=⋯=λs+N-=λ¯, corresponding eigenvectors are denoted by φs+1,φs+2,…,φs+N-. In the case of N->1, the calculation for the derivatives of eigenvalues is not straight forward. Based on a mathematical perturbation technique, a change Δρe of one arbitrarily selected design parameter ρe is firstly considered, and corresponding stiffness and mass matrices will be changed; i.e., (10)K+∂Κ∂ρeΔρeM+∂M∂ρeΔρe,e=1,…,NEThen the multiple eigenvalues and corresponding eigenvectors for the perturbed design can be written as(11)λ~jρe+εΔρe=λj+εujρe,Δρe+oεφ~jρe+εΔρe=φj+εvjρe,Δρe+oε,e=1,…,NEwhere uj and vj are unknown sensitivities of jth eigenvalue and corresponding eigenvector, respectively. Substituting (11) into (3), the first approximation can be obtained:(12)∂K∂ρe-λj∂M∂ρeφj+K-λjMvj=ujMφjNote that any linear combination of eigenvectors φ~i=∑k=ss+N-βkiφk,i=s,…,s+N- will satisfy (3), where βki are unknown coefficients. Premultiplying (12) by φ~iT has(13)φ~iT∂K∂ρe-λj∂M∂ρeφj=ujφ~iTMφj,j=s+1,…,s+N-and we have(14)∑k=ss+N-βkiφkT∂K∂ρe-λj∂M∂ρeφj-ujδij=0,i,j=s+1,…,s+N-,e=1,…,,NEFinally, the sensitivities uj=Δλj,j=s…,s+N- of multiple eigenvalues are determined by(15)detφkT∂K∂ρe-λj∂M∂ρeφj-Δλjδij=0,i,j=s+1,…,s+N-,e=1,…,,NFor 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 ∂λj/∂ρe as(16)∂λ~j∂ρe=1ρe∑i=1NEWl∑i=1NEWlρl∂λj∂ρeand the weight factor Wl in (16) is defined as(17)Wl=max0,rmin-diste,l,l∈NE∣diste,l≤rmin,e=1,…,NEwhere the operator dist(e,l) is defined as the distance between the center of element e and the center of element l, and rmin is the prescribed filter radius.


## 2.1. Topology Optimization Problems

The topology optimization problem for maximization of theith eigenfrequency of vibrating structures can be stated as [2] (1)maxρ1,…,ρNEβsubject to(2)β-ωl2≤0,l=i,i+1,…,J(3)Kφj=ωj2Mφj,j=1,…,J,(4)φjTMφk=δjk,j,k=1,…,J,,(5)∑e=1NEρeve≤V-,,(6)0≤ρ_≤ρe≤1,e=1,…,NE..where β is the lower bound for the ith eigenfrequencies. K and M are n×n the stiffness matrix and mass matrix, respectively. The ith eigenfrequency ωi and the corresponding eigenvector φi can be obtained by solving (3). J represents the number of the candidate eigenpairs, and it is suggested to be larger than i that prevents the ith eigenfrequency exchange its order with the neighbouring eigenfrequency during the optimization process [2]. Let these eigenfrequencies be ordered so that 0<ω1≤ω2≤⋯≤ωJ and the corresponding eigenvectors [φ1,φ2,…,φJ] satisfy constraint equation (4), where δjk is Kronecker’s delta. In (5) and (6), the design variable ρe(e=1,…,NE) is an element density and ρ_≤0.001 is used to prevent the mass and stiffness matrices from being singular. NE and ve denote the number of finite elements and the volume of element e, respectively. V¯ 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 paper(7)ke=15ρe5+ρe16ke0me=ρeme0where ke and me are the stiffness matrix and mass matrix of element e. ke0 and me0 represent element stiffness matrix and mass matrix of solid element e.


## 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ρe(e=1,…,NE) 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 thejth eigenfrequency ωj(λj=ωj2) is simple; i.e., ωj-1<ωj<ωj+1, and then the corresponding eigenvector φj is unique. By directly differentiating the constraint equation (3) with respect to ρe, we have(8)K-λjM∂φj∂ρe+∂K∂ρe-λj∂M∂ρe-∂λj∂ρeMφj=0,e=1,…NEPremultiplying (8) by φjT and combining with (3) and (4), the sensitivity of a simple eigenvalue λj can be obtained by(9)∂λj∂ρe=φjT∂K∂ρe-λj∂M∂ρeφj,e=1,…,NEwhere the derivatives of the matrices K and M can be obtained from (7).When the generalized eigenvalue problem in (3) has an N--fold multiple eigenvalues λj=ωj2,j=s+1,s+2,…,s+N-, with s≥0 and s+N-≤J, that is to say, λs+1=λs+2=⋯=λs+N-=λ¯, corresponding eigenvectors are denoted by φs+1,φs+2,…,φs+N-. In the case of N->1, the calculation for the derivatives of eigenvalues is not straight forward. Based on a mathematical perturbation technique, a change Δρe of one arbitrarily selected design parameter ρe is firstly considered, and corresponding stiffness and mass matrices will be changed; i.e., (10)K+∂Κ∂ρeΔρeM+∂M∂ρeΔρe,e=1,…,NEThen the multiple eigenvalues and corresponding eigenvectors for the perturbed design can be written as(11)λ~jρe+εΔρe=λj+εujρe,Δρe+oεφ~jρe+εΔρe=φj+εvjρe,Δρe+oε,e=1,…,NEwhere uj and vj are unknown sensitivities of jth eigenvalue and corresponding eigenvector, respectively. Substituting (11) into (3), the first approximation can be obtained:(12)∂K∂ρe-λj∂M∂ρeφj+K-λjMvj=ujMφjNote that any linear combination of eigenvectors φ~i=∑k=ss+N-βkiφk,i=s,…,s+N- will satisfy (3), where βki are unknown coefficients. Premultiplying (12) by φ~iT has(13)φ~iT∂K∂ρe-λj∂M∂ρeφj=ujφ~iTMφj,j=s+1,…,s+N-and we have(14)∑k=ss+N-βkiφkT∂K∂ρe-λj∂M∂ρeφj-ujδij=0,i,j=s+1,…,s+N-,e=1,…,,NEFinally, the sensitivities uj=Δλj,j=s…,s+N- of multiple eigenvalues are determined by(15)detφkT∂K∂ρe-λj∂M∂ρeφj-Δλjδij=0,i,j=s+1,…,s+N-,e=1,…,,NFor 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 ∂λj/∂ρe as(16)∂λ~j∂ρe=1ρe∑i=1NEWl∑i=1NEWlρl∂λj∂ρeand the weight factor Wl in (16) is defined as(17)Wl=max0,rmin-diste,l,l∈NE∣diste,l≤rmin,e=1,…,NEwhere the operator dist(e,l) is defined as the distance between the center of element e and the center of element l, and rmin 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 procedure(18)KΦ=MΦΛwhere K and M are the n×n stiffness and mass matrices, respectively. The p lowest eigenvalues Λ=diag(ω12,ω22,…,ωp2) and corresponding Φ=φ1,φ2,…,φp can be calculated via subspace iteration with shifting [23] or shifted Lanczos algorithm [24]. In addition, a series of LDLT factorizations of the shifted stiffness matrices K-μjMj=1,…,m,μj≥0 is involved in the solution procedure; m is the number of these triangular factorizations.According to the BCAS, a modified eigenvalue problem is introduced as follows:(19)KMΦM=MMΦMΛMwhere KM=K+ΔK, MM=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 μl can be given as (20)K-μlM+ΔK-μlΔMΦM=MMΦMΛ-Mwhere Λ-M=ΛM-μlI. The two eigenvalue problems in (19) and (20) have the same eigenvectors. It can be derived from (20) that(21)ΦM=I+K-μlM-1ΔK-μlΔM-1K-μlM-1MMΦMΛ-Mand (21) can be expanded as(22)ΦM=I-B+B2-…K-μlM-1MMΦMΛ-Mwhere(23)B=K-μlM-1ΔK-μlΔMSubstituting Φ for ΦM and dropping Λ-M in (22), the p×s basis vectors RB=[R1,R2,…,Rs] in the BCAS method are evaluated by(24)R1=K-μlM-1MMΦ(25)Rj=-BRj-1,j=2…,sAs the decomposed form ofK-μlM is usually given from the initial design, the calculation of RB only contains forward and backward substitutions. The stiffness matrix KM and mass matrix MM are condensed as(26)KR=RBTKMRB,MR=RBTMMRBThe modified analysis equation (19) can be approximated by the following reduced eigenvalue problem(27)KRΘR=MRΘRΛRwhere ΛR and ΘR are the matrices of p eigenvalues and the corresponding eigenvectors. The matrix of eigenvalues and corresponding matrix of eigenvectors for the modified analysis equation (19) can be evaluated by(28)ΦM=RBΘR,ΛM=ΛRFor BCAS method, the approximate solutions can be high accuracy achieved with s=2. By using the information of the eigenvalues for the initial design, the shifted frequency ul in (20) can be chosen as(29)μl-λ1+λp2=min1≤j≤mμj-λ1+λp2When the number of required eigenpairs is small, we usually set μl=0. 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 ρk and ρpE is defined as(30)ρk-ρpE≤εnormand the Modal Assurance Criterion (MAC) value is defined by (31)MACρk,ρpE=ρkTρpE2ρkTρkρpETρpEwhere ρk=[ρk,1,…,ρk,NE], ρpE=[ρp,1E,…,ρp,NEE] are the design variable vectors at kth step and the one corresponding to the previous exact solution, respectively.The convergence criteria [25] of the optimization procedure are set to (32)fk+1-fkfk≤10-4where fk is the objective value of kth step. Meanwhile, MAC(ρA,ρE) is used to illustrate the accuracy of the result, where ρA and ρE 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.Figure 1
The flow chart of the proposed optimization procedure.


## 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,J must be chosen to be sufficiently large to include all members of a possible N--fold eigenfrequency, such as J=j+2 or J=j+3. 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 8m×4m, and its thickness is 1m. 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 E=1.0×1010pa, Poisson’s ratio ν=0.3, and the mass density ρ=1000kg/m3. A concentrated nonstructural masse M=16000kg is applied in the corners of the right edge.Figure 2
A 2D cantilever beam with a concentrated nonstructural mass.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 Figure3. 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.Table 1
Results of maximization the first natural frequency of the 2D cantilever beam.
Method Design iterations Total SI iterations Total BCAS iterations Objective value MAC value CPU times (sec.) Exact method 54 223 0 77.09  129.6 Proposed method 53 86 130 76.35 0.982 83.9Figure 3
The converged results for a 2D cantilever beam: (a) exact method; (b) proposed method.
(a) (b)Figure 4
Convergence curve for the objective function of a 2D cantilever beam: proposed method versus exact method.


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

The beam of dimension14m×2m is clamped on both ends as shown in Figure 5, and its thickness is 0.01m. It is divided into 280×40 four-node quadratic elements with 22878 DOFs. Young’s module of material is set to E=2.5×1010pa, Poisson’s ratio ν=0.3, mass density ρ=2.5×103kg/m3. A concentrated nonstructural mass M=125kg is placed at the center.Figure 5
A 2D rectangle clamped beam with concentrated mass.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 Figure6. 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.Table 2
Results of maximization the first natural frequency of the 2D rectangle clamped beam.
Method Design iterations Total SI iterations Total BCAS iterations Objective value MAC value CPU times (sec.) Exact method 51 218 0 126.91  449.6 Proposed method 36 93 47 123.56 0.975 297.9Figure 6
The converged results for a 2D rectangle clamped beam: (a) exact method; (b) proposed method.
(a) (b)Figure 7
Convergence curve for the objective function of a 2D rectangle clamped beam: proposed method versus exact method.


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

The structure has the domain of20m×20m×1m 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 E=1.0×1011pa, Poisson’s ratio is ν=0.3, and the mass density is ρ=7800kg/m3. An additional artificial mass M=3.12×105kg is applied at the center of the structure.Figure 8
A 3D clamped square plate with concentrate mass.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 Figures9(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 ω2=ω3=65.37rad/s. 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%.Table 3
Results of maximization the second natural frequency of the 3D clamped square plate.
Method Design iterations Total SI iterations Total BCAS iterations Objective value MAC value CPU times (sec.) Exact method 31 122 0 177.89  2428.5 Proposed method 28 63 43 175.95 0.996 1652.0Figure 9
Optimal designs for a 3D clamped square plate: (a) exact method; (b) proposed method.
(a) (b)Figure 10
Maximizing the second eigenfrequency of a 3D clamped square plate. (a) Iteration histories of the first three eigenfrequencies obtained by exact method. (b) The comparison of convergence curve for the objective function based on exact method and the proposed method.
(a) (b)


## 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 8m×4m, and its thickness is 1m. 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 E=1.0×1010pa, Poisson’s ratio ν=0.3, and the mass density ρ=1000kg/m3. A concentrated nonstructural masse M=16000kg is applied in the corners of the right edge.Figure 2
A 2D cantilever beam with a concentrated nonstructural mass.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 Figure3. 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.Table 1
Results of maximization the first natural frequency of the 2D cantilever beam.
Method Design iterations Total SI iterations Total BCAS iterations Objective value MAC value CPU times (sec.) Exact method 54 223 0 77.09  129.6 Proposed method 53 86 130 76.35 0.982 83.9Figure 3
The converged results for a 2D cantilever beam: (a) exact method; (b) proposed method.
(a) (b)Figure 4
Convergence curve for the objective function of a 2D cantilever beam: proposed method versus exact method.


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

The beam of dimension14m×2m is clamped on both ends as shown in Figure 5, and its thickness is 0.01m. It is divided into 280×40 four-node quadratic elements with 22878 DOFs. Young’s module of material is set to E=2.5×1010pa, Poisson’s ratio ν=0.3, mass density ρ=2.5×103kg/m3. A concentrated nonstructural mass M=125kg is placed at the center.Figure 5
A 2D rectangle clamped beam with concentrated mass.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 Figure6. 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.Table 2
Results of maximization the first natural frequency of the 2D rectangle clamped beam.
Method Design iterations Total SI iterations Total BCAS iterations Objective value MAC value CPU times (sec.) Exact method 51 218 0 126.91  449.6 Proposed method 36 93 47 123.56 0.975 297.9Figure 6
The converged results for a 2D rectangle clamped beam: (a) exact method; (b) proposed method.
(a) (b)Figure 7
Convergence curve for the objective function of a 2D rectangle clamped beam: proposed method versus exact method.


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

The structure has the domain of20m×20m×1m 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 E=1.0×1011pa, Poisson’s ratio is ν=0.3, and the mass density is ρ=7800kg/m3. An additional artificial mass M=3.12×105kg is applied at the center of the structure.Figure 8
A 3D clamped square plate with concentrate mass.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 Figures9(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 ω2=ω3=65.37rad/s. 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%.Table 3
Results of maximization the second natural frequency of the 3D clamped square plate.
Method Design iterations Total SI iterations Total BCAS iterations Objective value MAC value CPU times (sec.) Exact method 31 122 0 177.89  2428.5 Proposed method 28 63 43 175.95 0.996 1652.0Figure 9
Optimal designs for a 3D clamped square plate: (a) exact method; (b) proposed method.
(a) (b)Figure 10
Maximizing the second eigenfrequency of a 3D clamped square plate. (a) Iteration histories of the first three eigenfrequencies obtained by exact method. (b) The comparison of convergence curve for the objective function based on exact method and the proposed method.
(a) (b)


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



---

*Source: 2728408-2018-11-05.xml*