A MAC based excitation frequency increasing method for structural topology optimization under harmonic excitations

This work is focused on the topology optimization of structures that are subjected to harmonic force excitation with prescribed frequency and amplitude. As an important objective of such a design problem, the natural resonance frequency of the structure is driven far away from the prescribed excitation frequency for the purpose of avoiding resonance and reducing the vibration level. Therefore when the excitation frequency is higher than the natural resonance frequency of the structure, the natural resonance frequency will decrease, then the optimum topology configuration will be distorted with large amount of gray elements. A MAC (Modal Assurance Criteria) based excitation frequency increasing method is proposed to obtain a desired configuration. MAC is adopted here to track the natural resonance frequency which can provide the baseline reference for the current excitation frequency during the optimum iterative procedure. Then the excitation frequency increases progressively up to its originally prescribed value. By means of numerical examples, the proposed formulation can generate effective topology configurations which can avoid resonance.


Introduction
Topology optimization has been recognized as an effective approach to figure out the structure layout during the conceptual design phase since the original idea of homogenization-based design method was proposed [1].Up to now, topology optimization has received remarkable success in both theoretical studies and practical applications [2][3][4], where topology optimization of continuum structures under dynamic loading is one challenge.Correlative researches in the field were mainly focused on two types of problems: one is topology optimization related to dynamic characteristics, the other is related to dynamic responses.
Natural frequencies and mode shapes are two major optimization objects in topology optimization related to dynamic characteristics.It follows some general objectives: The first objective is maximization of the specified structural eigenvalues such as fundamental or high-order natural frequencies [5]; secondly the gaps between the specified eigenvalues and frequencies must be maximized [6].The final objective is to optimize a structure to obtain prescribed desired natural frequencies and mode shapes [7].
Another type of problems is to minimize the dynamic responses such as displacement, velocity and acceleration amplitude or dynamic compliance in topology optimization related to dynamic responses.This paper mainly deals with the second case where the displacement amplitude under harmonic loads is involved.Physically, minimization of dynamic responses under harmonic loads is of great importance since harmonic vibration may be easily caused by rotating or reciprocating components.Besides, complicated periodic excitations can be transformed into a set of harmonic excitations by use of the Fourier series.
Relevant researches have already be carried out about topology optimization related to dynamic responses.Ma et al. [8] minimized structure dynamic compliance under harmonic loads using the homogenization method.Shu et al. [9] studied the topology optimization for minimizing frequency response using level set method.Yoon [10] conducted a comprehensive investigation for three model reduction schemes namely mode displacement method (MDM), Ritz vector method and quasi-static Ritz vector method in topology optimization under harmonic loads.
Xiang et al. [11] introduced a differentiable q-norm form of response function's peak value in dynamics response in truss sizing and shape optimization under bandwidth frequency excitation.Yang and Li [12] adopted mode tracking technique to minimize structural dynamic compliance at resonance frequencies in thermal environments.Liu et al. [13] made a comparative study about the accuracy and efficiency of harmonic responses calculating using mode displacement method, mode acceleration method and full method.
Although many considerable efforts have been made about the structure topology optimization under harmonic loads, the problem of distorted configuration with high-frequency excitation has not been well solved yet, as well as the boundary of low-frequency and high-frequency is not clear.Jog [14] presented that the topology optimization for minimizing displacement amplitude at a user-defined point in the structure may result in distorted configurations with large amount of intermediate density elements.Olhoff and Du [15] minimized the dynamic compliance of structure subject to harmonic loads, optimization results show that when the natural resonance frequency of the initial structure is less than the given excitation frequency, the natural resonance frequency will decrease and the static compliance of the structure will increase very quickly.Then he started out the design with a small value of excitation frequency and sequentially increased the value up to its prescribed value.It works well when the natural resonance frequency of the structure is always the same as the fundamental eigenfrequency.But unfortunately, the natural resonance frequency of most structures are not the same, and mode switching may occur during the optimization procedure which makes it more difficult.
This work focuses on the topology optimization with the excitation frequency higher than the natural resonance frequency.In order to obtain a clear anti-resonance configuration, a MAC based excitation frequency increasing method is proposed.In this way, the limit value of the natural resonance frequency can be estimated to distinguish low-frequency and high-frequency.

Optimization problem 2.1 Harmonic response analysis
As is known, the dynamic equilibrium equation of a discretized n-DOF (degree of freedom) structure under harmonic loads can be written as: where M, C, K are the mass matrix, damping matrix and stiffness matrix respectively.u(t) denotes the displacement vector, P(t) represents the harmonic external loading vector which can be expressed as PðtÞ ¼ Pe jxpt ðj 2 ¼ À1Þ, P and x p denote the magnitude vector of harmonic force and the given excitation frequency.
To implement the mode superposition method, the eigenfrequencies and eigenvectors are firstly obtained by solving the corresponding equation of free vibration as follows: The ith circular eigenfrequency x i and eigenvector U i can be obtained by solving the free vibration system characteristic equation: The mode shape matrix U = [U 1 , U 2 ,. . .U n ] is normalized by mass matrix, so the following relations hold: where n i denotes the ith damping ratio, Rayleigh damping leads to the following relation between damping ratio and frequency: a and b are two Rayleigh damping factors.By coordinating transformation, where y(t) is the vector of generalized coordinates.By substituting equation ( 6) into equation (7) and premultiplying U T , we obtain a set of n uncoupled equations of motion.
The solution of the above equation can be expressed as: Thus, the displacement response can be obtained by substituting equation ( 8) into equation ( 6): As for the above displacement response expression, if all n modes are taken into account, the result would the exact solution.However, considering the computing efficiency, only lower l modes are used to calculate the displacement response approximately, which is called mode displacement method (MDM).
When it comes to large-scale problems, a small l will generate big truncation errors, which may lead to an undesirable result.While, the mode acceleration method (MAM) is able to reduce errors by including the effects of the truncated modes using a static analysis [16].
The solution of the uncoupled motion equations in equation ( 7) can be rewritten as The use of equation (11) in equation ( 10) results in It is noted that, the inverse of the stiffness matrix can be represented by using all eigenmodes [17]: Notice that in MDM, just lower l modes are used for all the three terms in equation (12).But in MAM, we can obtain the first part exactly by solving a corresponding static problem and using equation ( 13) to include all n modes.As for the second and third parts in equation ( 12), it can be expressed as equation ( 14) according to equation ( 11): Hence, the displacement response approximated by MAM can be obtained through substituting equation ( 13) into equation ( 12) and combining with equation ( 14): An equivalent form of equation ( 15) can be obtained by using equation (13) in equation ( 15): Comparing between equation (10) and equation ( 16), the second term of equation ( 16) related to MAM takes the effects of higher order modes into account to some degree, so that the computing accuracy is improved apparently.Besselink et al. [17] made comparisons between MDM and MAM, he found that MAM outperformed the MDM in all cases.Hence, we adopt the MAM to calculate the displacement response.

Mode tracking
The MAC (Modal Assurance Criteria) has been used intensively in experimental modal analysis to measure the relevance between tested modes and calculated modes.It is adopted here to track the resonance mode successfully even though Mode switching occurs during the optimization procedure.The definition of MAC is (see [18]): In equation (17), U ref denotes the reference eigenvector which is the resonance mode shape of the initial structure; U cur,i refers to the ith eigenvector of the structure in current iteration step.The value of MAC varies from 0 to 1, when the two eigenvectors are orthogonal to each other, it equals 0; when it equals 1, the two vectors represent exactly the same mode shape.Generally the two modes cannot be the same, and the mode with the highest value of MAC is identified as the resonance mode shape, and the corresponding frequency is the natural resonance frequency in current iteration step.

MAC based structural topology optimization
under harmonic loads

Topology optimization formulation
The topology optimization problem for minimizing the displacement amplitude of concerned DOFs can be stated as: where g h is pseudo-density variable of element h, N d denotes the number of the design elements; g 0 is the lower bound of the pseudo-density variables, which is used to prevent the mass, stiffness and damping matrices from becoming singular, here g 0 = 0.001; ||u s (t)|| denotes the displacement amplitude of concerned DOFs.V and V U are the solid volume fraction and its upper bound, respectively.
It is known that in topological optimization of continuum structures subjected to dynamic or inertial loads, the interpolation scheme of Solid Isotropic Microstructure with Penalization (SIMP) would lead to localized modes because of the mismatch between element stiffness and mass when pseudo-density variable g i takes a small value.In fact, a variety of interpolation schemes have been proposed to eliminate the localized modes in some degree, such as Rational Approximation of Material Properties (RAMP) (Stolpe and Svanberg [19]), and Polynomial Interpolation Scheme (PIS) (Zhu et al. [20]).PIS is adopted here.
where, m h and k h denote the mass matrix and stiffness matrix of element h respectively.m h0 and k h0 are the mass and stiffness matrix of a corresponding element h with fully solid material.

Sensitivity analysis
Substitution of equation ( 8) into equation (15) gives Suppose Thus, the displacements can be express as The sensitivity of displacement response can be written as where, P represents the design independent external harmonic loads, K À1 P is the static displacement of the structure, which can be obtain by implementing an additional static analysis, its derivatives can obtained by means of adjoint method [21].As for the second term, it requires the sensitivities of eigenfrequencies and eigenvectors.Differentiating the free vibration system characteristic equation in equation (3) with respect to the design variable g h and premultiplying U T i yields Since the mass matrix and stiffness matrix are symmetric, the second part equals zero according to equation (3).Thus, we can obtain the sensitivity of eigenfrequency: The derivatives of the eigenvectors can be written as follows [22] where b ir is calculated as follows: The displacement amplitude of concerned DOFs in structure with damping yields where X s denotes the complex displacement.Then, we can derive the sensitivity of displacement amplitude through the chain rule:

Distorted configuration
A 3D cantilever beam is designed to illustrate the phenomenon of distorted configuration under high excitation frequency.The structure has a size of 0.8 m • 0.4 m • 0.04 mm which is clamped at the left side, as shown in Figure 1.The design domain is meshed with 80 • 40 • 4 solid elements.A harmonic force is applied at the middle node of the edge in the lower right corner with the amplitude of 5000 N. The parameters of material properties are as follows: E 0 = 2.1 • 10 11 Pa, l = 0.3, q 0 = 7800 Kg/m 3 , Rayleigh damping factors a = 10 À2 b = 10 À5 .The optimization objective is to minimize the vertical displacement amplitude at the loading position and the volume fraction is constrained to be less than 50%.
At first, we set the initial value of all variables as 0.5, then modal analysis is implemented to obtain the resonant frequency and reference eigenvector.The first four mode shapes and its national frequencies are shown in Figure 2, where f ni is the ith national frequency.As for the first three  mode shapes, they vibrate out of XY plane so that those natural mode will not be a resonance one.On the contrary, the fourth mode shape swings along the vertical direction, which makes most contributions to the target vertical displacement amplitude.As a conclusion, the fourth mode shape is the first resonance mode, corresponding national frequency f n4 is the initial resonance frequency f r0 , which can be also identify through frequency response analysis (Figure 3).Two harmonic forces with different level of given excitation frequency namely f p = 150 Hz < f r0 and f p = 160 Hz > f r0 are applied on the structure separately.The two optimized structures and their iteration history are presented in Figure 4.As we can see in Figure 4c, when the excitation frequency is less than the natural resonance frequency, minimization of the displacement amplitude drives the natural frequency away from the resonance point by increasing low national frequencies especially the natural resonance frequency, so that its static displacement ||u s ||(f P = 0) decreases, which means that the structure becomes stiffer.On the contrary, high excitation frequency results in decreasing of low national frequencies especially the natural resonance frequency to drive the natural resonance frequency away from the resonance point (Figure 4d).As a consequence, the optimization objective decreases a lot, but the static displacement increases, which means a weaker structure is obtained.It is known that a typical stiffness design is to find the optimal lay-out of a structure with the purpose of maximizing the performance of limit material.So that all materials distribute among the loading path, finally a clear and strong structure configuration is obtained such as Figure 4a, since its stiffness increases also.In contrast, because of the decrease of stiffness under high excitation frequency, optimization cannot motivate the material to distribute among the loading path, therefore it failed to generate a clear configuration (Figure 4b).Similar distorted configurations can be found in references [14,15].

MAC based excitation frequency increasing method
As is mentioned above, in order to obtain a clear configuration, it is critical to avoid the decrease of stiffness during optimization procedure.Because of the relation between the movement of natural resonance frequency and excitation   frequency, it is important to guarantee the value of excitation frequency below the natural resonance frequency.The implementation of MAC based excitation frequency increasing method can be stated as follows: Step 1, harmonic response analysis and modal analysis are implemented to determine the resonance frequency and extract the resonance mode as reference eigenvector with initial structure.
Step 2, MAC of each mode shape is evaluated to identify the resonance mode with the highest value of MAC, and the corresponding frequency is the natural resonance frequency f ri in ith iteration step.
Step 3, the excitation frequency f ei in ith iteration step can be determined according to the value of f ri : if the natural resonance frequency is equal or lesser than the prescribed frequency plus a little increment namely f ri f p + Df (f 0 ri in Figure 5), then f ei = f ri À Df; otherwise f ri > f p + Df (f 00 ri in Figure 5), then f ei = f p .Finally, f ei increases along with the raising of the natural resonance frequency f ri to the prescribed value.The completed optimization strategy is represented in Figure 6.
Since the excitation frequency is always lower than the natural resonance frequency during the optimization procedure, the final optimized configuration is fine without distortion.

Numerical examples
The same 3D cantilever beam (Figure 1) under excitation frequency f p = 480 Hz > f r0 is designed to illustrate the validity of the proposed design strategy in this section.As shown in Figure 7, a clear optimized configuration is obtained after 128 iterations.The optimization procedure can       Figure 9 shows the iteration history of the resonance mode order and the max value of MAC.As we can see, the resonance mode order changes during the optimization procedure, but it is tracked successfully by identifying the max value of MAC.The max value of MAC holds a value above 0.9, which indicates that large deformation will not occur on the resonance mode shape.
For the purpose of comparison, the same structure is designed with different optimization objectives.The first one is the integral of displacement amplitude in a frequency interval (f A , f B ) written as equation (30).A clear optimized configuration (Figure 10a) is obtained by means of the same integral calculation method in reference [13].The second is the static compliance under same displacement amplitude force, Figure 10b shows the results.Compared to the initial structure, the static displacements of three optimized structures all decrease dramatically (Table 1), which demonstrate the significance of structural stiffness in attaining a clear configuration.
Min : Figure 11 shows the displacement response curves of three different configuration.As we can see, the optimized configuration in static condition is failed to avoid resonance that a major peak value will appear in the frequency interval     (0, 480 Hz).As for the optimized configuration in a frequency interval, the value of the optimization objectives ||u s (t)|| and R 480 0 jju s jjdf are both the minimum one in all the three optimized configurations (Table 1).It is mainly because that the peak values of displacement decrease and troughs appear in the frequency interval.Even though, resonance will occur when external excitation reach the exact value.Luckily the optimized configurations using the proposed method can avoid resonance successfully by increasing the natural resonance frequency larger than the prescribed value of excitation frequency.
It should be noted that the natural resonance frequency has a limit value, when the prescribed excitation frequency exceeds this limit, the proposed method is inapplicable.To detect this limit value, we keep the excitation frequency always lower than the natural resonance frequency in the optimization procedure.Figure 12 shows the iteration history of the natural resonance frequency, it is observed that when the excitation frequency is closed to the limit value, optimization objective of displacement amplitude can decrease only by reducing the resonance frequency, since it is the only one path to drive the excitation frequency away from the natural resonance frequency.Then the natural resonance frequency will increase again along with the increment of the excitation frequency.Finally, the resonance frequency fluctuates around the limit value.Therefore we can obtain this limit value approximately.This significant value can be used as the boundary of lowfrequency and high-frequency excitation problems.In this way, the range of low-frequency extends a lot compared with the traditional criterion which just simply sets the natural resonance frequency of initial structure as the boundary line.It means that we can obtain clear optimum anti-resonance structures in a greater scope by employing the proposed method.As for the high-frequency excitation problems, its optimization objective may need a change or additional constraints should be taken into account to obtain desired results.

Conclusions
Topology optimization of structures with damping that are subjected to harmonic force excitation with prescribed frequency higher than its natural resonance frequency is carried out in this paper.A MAC based excitation frequency increasing method is proposed to avoid distorted configuration.Compared with other two configurations which use the integral of displacement amplitude in a frequency interval or the same displacement amplitude under static condition as optimization objective, this configuration has a significant advantage of avoiding resonance since its natural frequency increases above the given excitation frequency.In addition, the limit value of the natural resonance frequency can be calculated approximately to distinguish low-frequency and high-frequency excitation problems so that the range of low-frequency extends a lot.As a result, we can obtain clear optimum anti-resonant structures in a greater scope by employing the proposed method.

Figure 1 .
Figure 1.3D cantilever beam subjected to a harmonic force.

6 T.
Liu et al.:  Int.J. Simul.Multisci.Des.Optim.2017, 8, A4 be divided into two stages as shown in Figure8: in the first stage, the excitation frequency increases rapidly along with the natural resonance frequency to the prescribed value, the displacement amplitude moves irregularly due to the change of excitation frequency.A small frequency increment Df accelerate the increase of natural resonance frequency, which means that raising the value of natural resonance frequency to the prescribed one needs fewer iteration steps.Moreover it will not affect the final optimization results.Here Df is set at the value of 15 Hz.In the second one, the excitation frequency remain the prescribed value unchanged, after a small fluctuation, the natural resonance frequency increases smoothly, the displacement amplitude decrease smoothly also.

Figure 8 .
Figure 8. Iteration history of the displacement amplitude the excitation frequency along with the first resonance frequency.

Figure 9 .Figure 6 .
Figure 9. Iteration history of the resonance mode order and the max value of MAC.

Figure 10 .
Figure 10.Optimized configurations with different optimization objectives.(a) Optimized configuration in a frequency interval.(b) Optimized configuration in static condition.

Figure 11 .
Figure 11.Comparison of displacement amplitudes of three optimization structures.

Figure 12 .
Figure 12.Iteration history of the resonance frequency.

Table 1 .
The displacement values of different structures.