© T. Liu et al., published by EDP Sciences, 2017

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

2 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:(1)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 and ω _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:(2)

The ith circular eigenfrequency ω _i and eigenvector Φ _i can be obtained by solving the free vibration system characteristic equation:(3)

The mode shape matrix Φ = [Φ ₁, Φ ₂,… Φ _n ] is normalized by mass matrix, so the following relations hold:(4)where ξ _i denotes the ith damping ratio, Rayleigh damping leads to the following relation between damping ratio and frequency:(5) α and β are two Rayleigh damping factors.

By coordinating transformation,(6)where y(t) is the vector of generalized coordinates. By substituting equation (6) into equation (7) and premultiplying Φ ^T , we obtain a set of n uncoupled equations of motion.(7)

The solution of the above equation can be expressed as:(8)

Thus, the displacement response can be obtained by substituting equation (8) into equation (6):(9)

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).(10)

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(11)

The use of equation (11) in equation (10) results in(12)

It is noted that, the inverse of the stiffness matrix can be represented by using all eigenmodes [17]:(13)

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):(14)

Hence, the displacement response approximated by MAM can be obtained through substituting equation (13) into equation (12) and combining with equation (14):(15)

An equivalent form of equation (15) can be obtained by using equation (13) in equation (15):(16)

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.

2.2 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]):(17)

In equation (17), Φ _ref denotes the reference eigenvector which is the resonance mode shape of the initial structure; Φ _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.

3 MAC based structural topology optimization under harmonic loads

3.1 Topology optimization formulation

The topology optimization problem for minimizing the displacement amplitude of concerned DOFs can be stated as:(18)where η _h is pseudo-density variable of element h, N _d denotes the number of the design elements; η ₀ is the lower bound of the pseudo-density variables, which is used to prevent the mass, stiffness and damping matrices from becoming singular, here η ₀ = 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 η _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.(19)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.

3.2 Sensitivity analysis

Substitution of equation (8) into equation (15) gives(20)

Suppose(21)

Thus, the displacements can be express as(22)

The sensitivity of displacement response can be written as(23)where, P represents the design independent external harmonic loads, K ⁻¹ 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 η _h and premultiplying yields(24)

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:(25)

The derivatives of the eigenvectors can be written as follows [22](26)where β _ir is calculated as follows:(27)

The displacement amplitude of concerned DOFs in structure with damping yields(28)where X _s denotes the complex displacement. Then, we can derive the sensitivity of displacement amplitude through the chain rule:(29)

3.3 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 ₀ = 2.1 × 10¹¹ Pa, μ = 0.3, ρ ₀ = 7800 Kg/m³, Rayleigh damping factors α = 10⁻² β = 10⁻⁵. 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%.

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

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

Figure 2. The first four mode shapes and natural frequencies of initial structure. (a) First mode shape, f n1 = 18.6 Hz. (b) Second mode shape, f n2 = 77.3 Hz. (c) Third mode shape, f n3 = 115.0 Hz. (d) Fourth mode shape, f n4 = 156.4 Hz.

Figure 3. Displacement amplitude of initial structure.

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

Figure 4. The optimized structures and iteration history under different excitation frequencies. (a) Optimized configuration (f p  = 150 Hz). (b) Optimized configuration (f p  = 160 Hz). (c) Iteration history (f p  = 150 Hz). (d) Iteration history (f p  = 160 Hz).

3.4 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  + Δf ( in Figure 5), then f _ei  = f _ri  − Δf; otherwise f _ri  > f _p  + Δf ( in Figure 5), then f _ei  = f _p .

  Figure 5. Excitation frequency depends on resonance frequency.

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

  Figure 6. Flow chart of the optimization strategy.

Since the excitation frequency is always lower than the natural resonance frequency during the optimization procedure, the final optimized configuration is fine without distortion.

4. 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 be divided into two stages as shown in Figure 8: 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 Δf 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 Δf 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 7. Optimized configuration with given excitation frequency (f p  = 480 Hz).

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

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.

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

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.(30)

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

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

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

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

Figure 12. Iteration history of the resonance frequency.

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

Acknowledgments

This work is supported by National Natural Science Foundation of China (11432011, 51521061, NSAF U1330124), the 111 Project (B07050), the Fundamental Research Funds for the Central Universities (3102014JC02020505).

References

All Tables

All Figures

	Figure 1. 3D cantilever beam subjected to a harmonic force.
In the text

	Figure 2. The first four mode shapes and natural frequencies of initial structure. (a) First mode shape, f n1 = 18.6 Hz. (b) Second mode shape, f n2 = 77.3 Hz. (c) Third mode shape, f n3 = 115.0 Hz. (d) Fourth mode shape, f n4 = 156.4 Hz.
In the text

	Figure 3. Displacement amplitude of initial structure.
In the text

	Figure 4. The optimized structures and iteration history under different excitation frequencies. (a) Optimized configuration (f p  = 150 Hz). (b) Optimized configuration (f p  = 160 Hz). (c) Iteration history (f p  = 150 Hz). (d) Iteration history (f p  = 160 Hz).
In the text

	Figure 5. Excitation frequency depends on resonance frequency.
In the text

	Figure 6. Flow chart of the optimization strategy.
In the text

	Figure 7. Optimized configuration with given excitation frequency (f p  = 480 Hz).
In the text

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

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

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

	Figure 11. Comparison of displacement amplitudes of three optimization structures.
In the text

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