A finite element based optimal vibration suppression for constrained layer damped rotating plates

1. Introduction

Constrained layer damping (CLD) is an effective approach for vibration suppressions. The fundamental mechanism of CLD is that the relative motion of two face layers leads to the deformation of the damping layer; the energy consumed through the deformation of the damping layer results in the damping effects. Relatively large damping ratios can be achieved through a wide range of available viscoelastic materials.

In literature, investigations have been conducted on constrained layer damping for beam structures. These researches widely employed an assumption proposed by Ross, Ungar and Kerwin (RKU) [1], suggesting that there is only shear deformation (pure shear) in the damping layer. In this assumption, the two face layers are treated as Euler-Beams while the core layer only experiences shear deformation. The motions associated with this assumption are illustrated in Fig. 1. For the convenience of the discussion, two face layers are called the constraining layer and the base layer respectively in this paper, while the middle layer is called the core layer. The layers are also denoted by numbers 1, 2 and 3, respectively, throughout this work. The displacements $u_1$ and $u_3$ are employed to represent the longitudinal displacements of face layers 1 and 3 respectively. $w$ is utilized to denote the transverse displacement of sandwich structure, which is identical for each of the three layers. Shear deformation in the core layer is determined from the relative motion between $u_1$ and$u_3$.

While the assumptions described above are widely used, they are not always valid in certain structures. Rao [2] studied the problem using the same transverse displacement assumption as that in [1] but with additional assumptions that the longitudinal displacement varies linearly across the thickness of each layer, and the displacement is continuous across the layer connections. In that work, all of three layers are treated as Timoshenko beams instead of Euler beams used in the RKU model. Zapfe and Lesieutre [3] presented a beam finite element that can be regarded as the finite element implementation of the Rao’s model. Based on the work of [3], Xie and Shepard [4] developed a new plate finite element to consider the compression effect of CLD structures.

Fig. 1RKU model displacements for a CLD sandwich beam

Recently the constrained layer damping is employed for the vibration suppressions of rotating structures, which has wide applications in industries such as wind turbines, helicopters, and vessels. Rotating beams, a simplification of rotating plates, have been extensively investigated, including fundamental dynamics [5] and constrained layer damping studies [6-7]. Fung et al. [8] studied active constrained layer damping (ACLD) application on rotating beam. Yoo et al. [9-11] studied the dynamics of rotating bare plates using analytical approach. Based on the work in [8], Liu et al. [12] developed a plate finite element model for a rotating bare plate. Liu et al. [13] extended their own work in [12] to the active control of rotating bare plates through the finite element approach. It is worth noting that all of the work mentioned above utilized the same assumptions as those in RKU model that there is only shear strain in the core layer, whereas the constraining and base layers experience only longitudinal extension strain.

Since it is often desired to obtain the best performance for a particular application, design optimization can often be a significant issue. Many optimization methods are available, such as calculus-based methods like the conjugate gradient method and differential correction. In some cases, the model associated with optimization issue is complicated or even unknown; it is difficult for these methods to handle such optimization issues. Genetic algorithm (GA) could be a good choice because GA is an optimization method based on stochastic search, requiring only the knowledge of the objective function and not any of its derivatives. GA uses natural evolutionary principles to obtain offspring with characteristics better than those of the parents. The optimization results are typically at or close to the global optimum. Some GA optimization methods have been applied to the constrained-layer damped structures. Hau et al. [14] constructed a FEM model which contains the frequency-dependent viscoelastic material for a rotating beam. This model was then combined with GA to obtain an optimum configuration.

In this research, a plate finite element model is originally developed to extend the previous research from the beam to the plate and is employed to model the multi-layer rotating plates. A GA approach is then used with the model to find the optimal thickness for both the damping layer and the constraining layer for vibration suppressions. A comparison is provided in the numerical examples, demonstrating that the optimum solutions obtained from GA can obtain the best vibration suppression effects while minimizing the utilization of the damping material. The proposed plate finite element model is applicable for multi-layer (more than three) rotating plates, and the design approach can be extended to optimize more complex damped rotating structures such as wind turbines in the future.

2. Finite element modeling and validation

2.1. Assumptions and kinematics relations

Prior to introducing the rotating plate finite element model, the assumptions as those in [2] should be introduced. Note that the assumptions in [2] are for beams.

1) The transverse displacements of all three layers are equal.

2) The longitudinal displacement is linearly distributed across the thickness of each layer.

3) There is no slip between the layers.

These assumptions are used in this work to model the rotating plate. Fig. 2 shows the proposed displacement field for the rotating plate with angular velocity $\omega$. $u_i$ and $v_i$ are the longitudinal displacements at the contacting interfaces of different layers, and the longitudinal displacement of each layer varies linearly between the top and the bottom surfaces as shown in Fig. 2. $w$ is the transverse displacement for all of the three layers.

Fig. 2Proposed displacement field

a)

b)

c)

Note that there are two coordinate systems in Fig. 2, i.e., a rotating coordinate system $oxyz$ attached to the rotating plate, and a global coordinate system $OXYZ$ that is fixed without rotation. The plate rotates with the rotating coordinate system $oxyz$ around the $Y$ axis of the fixed coordinate system $OXYZ\text{.}$ The transformation matrix between two coordinate systems is defined as [10]:

1

$A=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\theta & 0& -\mathrm{s}\mathrm{i}\mathrm{n}\theta \\ 0& 1& 0\\ \mathrm{s}\mathrm{i}\mathrm{n}\theta & 0& \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right],$

where $\theta$ is the angle between the $x$ axis and $X$ axis. Taking an arbitrary point on the rotating plate which moves from $P$ to $P\text{'}$ after the deformation, the:

$OP{\text{'}}_{xyz}=\left[\begin{array}{c}x\\ y\\ 0\end{array}\right]+\left[\begin{array}{c}u\\ v\\ w\end{array}\right]+\left[\begin{array}{c}-{\int }_0^x{\left(\frac{\partial w}{\partial x}\right)}^{2}dx\\ -{\int }_0^y{\left(\frac{\partial w}{\partial y}\right)}^{2}dy\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\right].$

The term: $\left[\begin{array}{c}-{\int }_0^x{\left(\frac{\partial w}{\partial x}\right)}^{2}dx\\ -{\int }_0^y{\left(\frac{\partial w}{\partial y}\right)}^{2}dy\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\right]$, is the so-called ‘stiffening effect’ [11], which are the coupling terms between the transverse displacement and the longitudinal displacement due to the centrifugal force of the rotation. The $u$, $v$, and $w$ are the linear displacement of point $P$ in the rotating coordinate system $oxyz$. The transformation of the displacement from the $oxy$z coordinate system to the $OXYZ$ coordinate system can be made as follows:

2

$O{P}_{XYZ}^{\text{'}}=A\left(\left[\begin{array}{c}x\\ y\\ 0\end{array}\right]+\left[\begin{array}{c}u\\ v\\ w\end{array}\right]+\left[\begin{array}{c}-{\int }_0^x{\left(\frac{\partial w}{\partial x}\right)}^{2}dx\\ -{\int }_0^y{\left(\frac{\partial w}{\partial y}\right)}^{2}dy\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\right]\right).$

With the proposed displacement field, the kinetic and strain energy can be obtained respectively. The equations of motion can then be derived accordingly, and they are presented in the following sections.

2.3. Finite element discretization

Once the kinetic and strain energy are obtained, the displacement field can be discretized using the new plate finite element. Fig. 3 shows the $i$th layer of the proposed new plate finite element. It is worth noting that the middle nodes are employed to avoid shear locking [3]. As a consequence, there are 25 degrees of freedom for the $i$th layer as shown in Fig. 3 and 41 degrees of freedom for all of three layers.

Fig. 3Proposed new plate finite element

a)

b) For nodes p, q, m, k

c) For all other nodes

The displacement vector of the $i$th layer with this new plate finite element becomes as follows:

6

$q_i={\left[\begin{array}{l}{u}_{pi}{u}_{p\left(i+1\right)}{v}_{pi}{v}_{p\left(i+1\right)}{w}_p{w}_1{u}_{qi}{u}_{q\left(i+1\right)}{v}_{qi}{v}_{q\left(i+1\right)}{w}_q{w}_2{w}_3{w}_4\\ u_{mi}{u}_{m\left(i+1\right)}{v}_{mi}{v}_{m\left(i+1\right)}{w}_m{w}_5{u}_{ki}{u}_{ki}{v}_{ki}{v}_{k\left(i+1\right)}{w}_k\end{array}\right]}^T.$

For the $i$th layer of plate finite element, the proposed displacement field (Fig. 3) can be organized as:

7

${\left[\begin{array}{c}u\\ v\\ w\end{array}\right]}_i=B\left[\begin{array}{c}{u}_i\\ u_{i+1}\\ v_i\\ v_{i+1}\\ w\end{array}\right]=B\left[\begin{array}{c}{F}_1\\ F_2\\ F_3\\ F_4\\ F_5\end{array}\right]q_i=BF{q}_i=N{q}_i,a,$

and the dimensions of matrix $B$ and $F$ are $\text{3×5}$ and $\text{5×25}$, respectively. Both matrixes are given as follows:

8

$B=\left[\begin{array}{ccccc}\frac{z}{H_i}& 1-\frac{z}{H_i}& 0& 0& 0\\ 0& 0& \frac{z}{H_i}& 1-\frac{z}{H_i}& 0\\ 0& 0& 0& 0& 1\end{array}\right],$

where $z$ and $H_i$ are respectively the transverse coordinate and the thickness of the $i$th layer. And:

9

$F=\left[\begin{array}{ccccccccccccccccccccccccc}{n}_1& 0& 0& 0& 0& 0& n_2& 0& 0& 0& 0& 0& 0& 0& n_3& 0& 0& 0& 0& 0& n_4& 0& 0& 0& 0\\ 0& n_1& 0& 0& 0& 0& 0& n_2& 0& 0& 0& 0& 0& 0& 0& n_3& 0& 0& 0& 0& 0& n_4& 0& 0& 0\\ 0& 0& n_1& 0& 0& 0& 0& 0& n_2& 0& 0& 0& 0& 0& 0& 0& n_3& 0& 0& 0& 0& 0& n_4& 0& 0\\ 0& 0& 0& n_1& 0& 0& 0& 0& 0& n_2& 0& 0& 0& 0& 0& 0& 0& n_3& 0& 0& 0& 0& 0& n_4& 0\\ 0& 0& 0& 0& b_1& b_2& 0& 0& 0& 0& b_3& b_4& b_5& b_6& 0& 0& 0& 0& b_7& b_8& 0& 0& 0& 0& b_9\end{array}\right],$

where $n_1~n_4$ are the shape functions for the longitudinal displacement and $b_1~b_9$ are the shape functions for the transverse displacement. Since they can be calculated according to many finite element textbooks, they are not given in detail in this work.

Also, the displacement vector $q_i$ can be assembled to the global displacement vector $q$ of the whole plate finite element with three layers:

10

$q=\sum _{i=1}^{3}T{R}_i{q}_i,$

where $T{R}_i$ is a $n\times \text{25}$ mapping matrix of the $i$th layer, and $n$ is the total number of degrees of freedom in the plate finite element (41 for three layers).

Substituting Eq. 7 into Eqs. (2), (3) and (5), and neglecting high order terms, the kinetic and strain energy of the $i$th layer can be reorganized as:

11

$V_i=\frac{\text{1}}{\text{6}}{\iiint }_{\begin{array}{c}{V}_i\end{array}}\left(E_i{q}_i^T\left({\left(\frac{\partial F_1}{\partial x}\right)}^T\left(\frac{\partial F_1}{\partial x}\right)+{\left(\frac{\partial F_1}{\partial x}\right)}^T\left(\frac{\partial F_1}{\partial x}\right)+{\left(\frac{\partial F_2}{\partial x}\right)}^T\left(\frac{\partial F_2}{\partial x}\right)\right)q_i+E_i{q}_i^T\left({\left(\frac{\partial F_3}{\partial y}\right)}^T\left(\frac{\partial F_3}{\partial y}\right)+{\left(\frac{\partial F_3}{\partial y}\right)}^T\left(\frac{\partial F_4}{\partial y}\right)+{\left(\frac{\partial F_4}{\partial y}\right)}^T\left(\frac{\partial F_4}{\partial y}\right)\right)q_i+\text{3}{G}_i{q}_i^T{\left(\frac{\partial F_5}{\partial x}-\frac{F_1-F_2}{\text{2}{H}_i}\right)}^T\left(\frac{\partial F_5}{\partial x}-\frac{F_1-F_2}{\text{2}{H}_i}\right)q_i+\text{3}{G}_i{q}_i^T{\left(\frac{\partial F_5}{\partial y}-\frac{F_3-F_4}{\text{2}{H}_i}\right)}^T\left(\frac{\partial F_5}{\partial y}-\frac{F_1-F_2}{\text{2}{H}_i}\right)q_i\right)dv,$

12

$T_i=\frac{1}{2}{\rho }_i{\iiint }_{\begin{array}{c}{V}_i\end{array}}\left({\dot{q}}_i^T{N}^{T}N{\dot{q}}_i+\omega \left(2x{N}_3{q}_i+{\dot{q}}_i^T{N}_3^T{N}_1{q}_i+q_i^T{N}_1^T{N}_3{\dot{q}}_i^T-{\dot{q}}_i^T{N}_1^T{N}_3{q}_i-{\dot{q}}_i^T{N}_3^T{N}_1{q}_i\right)+{\omega }^2\left(x^2+q_i^T{N}_1^T{N}_1{q}_i+q_i^T{N}_1^T{N}_1{q}_i\right)-x{q}_i^T{\int }_0^x{\left(\frac{\partial N_3}{\partial x}\right)}^T\left(\frac{\partial N_3}{\partial x}\right)dx{q}_i+2x{N}_1{q}_i\right)dv,$

where $N_i$ and $F_i$ are the $i$th row of the matrix $N$ and $F$, respectively.

3. Optimization via genetic algorithm (GA)

GA is not new to the field of mechanical designs. In 1975, Holland published the book of ‘Adaptation in Natural and Artificial System’ [15] and systematically described the topic of GA as applied to biological research. Since then, GA has been widely used in many fields such as structure optimizations. In general, GA uses stochastic search and natural evolutionary principles as the basis for optimizations. The algorithms first create a population of “individuals”, typically the binary strings coded to represent the parameters of the problem. The best elements of the population are then selected for various operations as shown in Fig. 4, and details on each operation can be found in [16]. The goal of each step in the algorithm is to obtain a new generation, referred to as offspring, which is closer to the optimal solution than the parent generation.

In this work, the thickness of the damping layer and the thickness of the constraining layer are treated as the design variables encoded into a binary string using a multi-parameter, mapped, fixed-point coding [16]. Each of the thicknesses is coded using 20 bits, so the total length of the string needed to describe the configuration of the sandwich beam is 40 bits. The population of each generation is composed of a specified number of such strings. The corresponding dynamic response of the rotating plate, associated with each string is used to calculate the fitness of the string. The fitness is used to ultimately determine the survivability of the particular string in the evolution. In this research, the fitness is defined as the summation of reciprocals of maximum amplitudes of the structural response at its first two resonant frequencies.

Fig. 4Procedure for genetic algorithm

The FEM code reads the strings generated from the GA code and calculates the dynamic response for each configuration. The dynamic response is then used in the embedded objective function in order to calculate the fitness for that particular solution. The particular objective functions considered for the optimization used herein are discussed below.

4. Numerical example

In order to examine the capabilities of the code, one can consider a constrained layer damped rotating plate with parameters as given in Table 2. For the GA parameters, the population size is 40 with single point crossover. Furthermore, forty bits are allocated for strings which represent the two design variables with twenty bits for each. The mutation rate is set as the reciprocal of the string length, which is 0.025. For this study, the first two mode amplitudes of driving point mobility at one selected point of the rotating plate will be minimized. The location of the selected point is at the point with one half of the plate length away from the free end and two third of the plate width away from the front side edge. Consequently, the response is presented in the frequency domain. While the driving point mobility amplitudes are minimized here, the method is flexible to handle other objective function.

Depending on the goals of the design, one can specify different objective functions for use in the GA optimization. Nakra [17] tried to maximize the system loss factor for certain modes. Zheng et al. [18] set the minimization of the vibration energy as the objective. Xie et al. [19] adopted GA to reduce the vibration amplitude of first two bending modes of composite beams. Since the response at the resonant frequencies is typically the most harmful, the initial objective function used here involves minimizing the driving point mobility amplitudes that would result from the simultaneous excitation of the first two modes (the first bending and the first torsion mode) at their respective resonance frequencies. For this case, the fitness (FS) is defined as $FS=\frac{1}{R_1}+\frac{1}{R_2}$ so that $R_1$ and $R_2$ can be minimized when the GA finds the configuration of the maximum fitness. Here, $R_1$ and $R_2$ are the first two mode amplitudes of the driving point mobility, respectively, at the selected point on the plate. As noted earlier, there are two design variables: 1) the thickness of damping layer $H_2$ and 2) the thickness of constraining layer $H_1$. In order to make the optimization practical, it is also necessary to specify the limits on the thickness of the two layers. For this problem, the thickness for these two layers must be kept within a predefined practical range, which can be expressed as:

where the maximum and minimum values of both layers are $H_{1\max }=H_{2\max }=$0.0012 m and $H_{1\min }=H_{2\min }=$0.0002 m, respectively. The decoded values of the binary strings are mapped into the range between the upper and lower limits automatically.

By using the above objective function and the layer thickness limits, the GA algorithm can find the optimal thicknesses for the two layers. In this case study, where the first two modes are considered simultaneously, and viewed as the criteria of the optimum configuration.

With the above settings, the GA and the FEM work together effectively such that the population appears to converge after the evolution of about sixteen generations, as shown in Fig. 5. The optimization takes about 72 hours (three days roughly) running on a 2.67 GHz PC, and most of the time is spent for the calculation of structure stiffness matrix which contains many differential operations. If a larger population size and more generations are used, the convergence of optimization process will take longer time.

Fig. 6 shows the responses of several different configurations. In the optimum configuration, the resulting thicknesses for the layers are 0.0012 m, 0.00085 m and 0.005 m for the base layer, the damping layer and the constraining layer, respectively. The responses of a bare rotating plate with a thickness of 0.005 m are used for the comparison purposes. The responses for the two configurations that respectively have the maximum and minimum layer thicknesses (i.e. $H_1=H_{1\mathrm{m}\mathrm{a}\mathrm{x}}$, $H_2=H_{2\mathrm{m}\mathrm{a}\mathrm{x}}$ and $H_1=H_{1\mathrm{m}\mathrm{i}\mathrm{n}}$, $H_2=H_{2\mathrm{m}\mathrm{i}\mathrm{n}}$) are also employed for more comparison. Fig. 6 shows that the optimum configuration performs about as well as the maximum thickness configuration and the responses for these two configurations are clearly much lower than the configuration without damping.

The reason that this work selects upper and lower end values of the thickness range is that engineers could benefit from this optimization by having maximum vibration reduction with the minimum material cost. It is not necessarily that more damping material can lead to better performance. As indicted previously, because of the complex modulus of viscoelastic material, the damping which to lower the vibration amplitude depends on the deformation of the damping layer. The more deformation in the damping layer, the more damping can be “squeezed” out of it. And this is the fundamental principle for the constrained layer damped structures. So, based on this consideration, following comparisons will made in order to illustrate the performance of proposed approach including CLD modeling, application, and optimization.

Fig. 5Evolution of the best-so-far fitness vs. generation, during the optimization process

Fig. 6Comparison on frequency response function (FRF) of different configurations: solid line – optimum configuration; dashed line – maximum thickness configuration; dash dot line – minimum thickness configuration; dotted line – bare rotating plate

Fig. 7 and 8 show detailed comparisons at the first and the second mode, respectively. It can be seen that the optimum configuration has about the same response level as the configurations with the maximum thickness at both resonant frequencies, while the damping layer thickness in the optimum configuration is lower than the configuration with the maximum thickness, indicating that less material can achieve about the same performance. Actually, the fitness values are 6597 for the optimum configuration, and 6471 for the maximum thickness configuration. Also, both configurations have a much lower response level than the minimum thickness configuration and the bare rotating plate at both resonant frequencies.

From the above comparisons, it can be seen that the first mode of the damped beam has a higher response than the response for the second mode. Since it is common to excite the first resonant frequency, the optimization of the first mode is more important and may be more important than the alternative selection of the optimum configuration over multiple resonant frequencies. Nevertheless, like some existing works [19], the selection of the best optimum approach must be based on some knowledge of the excitation and the critical response characteristics.

Fig. 7Comparison on frequency response function (FRF) at the first resonant frequency: solid line –optimum configuration; dashed line – maximum thickness configuration; dash dot line – minimum thickness configuration; dotted line – bare rotating plate

Fig. 8Comparison on frequency response function (FRF) at the second resonant frequency: solid line – optimum configuration; dashed line – maximum thickness configuration; dash dot line – minimum thickness configuration; dotted line – bare rotating plate

5. Conclusion

This research is the first attempt at using a plate finite element method to model the rotating plates with three layers. The proposed displacement field and degrees of freedom enable this finite element model to capture the shear and extension deformation in all three layers. This model is validated using the results of natural frequencies and damping ratios from the published literature. In order to find the geometric configurations to achieve the best desired structural response, the GA method is originally used to minimize the driving point mobility amplitudes of the constrained layer damped rotating plate in a frequency range that includes the first bending and the first torsion modes. The optimization result shows that GA can effectively work with the developed single layer plate finite element for rotating CLD plate to yield a satisfactory and very close vibration reduction to the configurations with a thicker layer of damping materials.