Locally resonant materials are widely used to improve the vibrational behavior of light-weight structures. Such materials can be created by integrating small beam-like structures into the cells of honeycomb structures. For the simulation of such materials a geometrically exact representation of the beam-like resonators is inefficient. Only the first one or two modes of the resonators influence the dynamic response of the host structure. Therefore it is beneficial, in order to reduce computation time, to replace the exact models of the resonators with SDOF or two-DOF systems which are capable of representing the first one or two modes of the resonators.

Honeycomb Structure

The following figure shows the resonator which is used for this project. It consists of a beam with a mass attached to its end. Because of limitations in the production the transition from the beam to the mass has to be inclined. The inclination helps further to avoid high stress peaks which would occur at sharp edges. 

beam-like resonator

The resonators are attached to the host structure which for this purpose is being divided into unit-cells of equal length. Each cell contains a resonator which is later modeled as SDOF or two-DOF system as already mentioned.

Unit cell with SDOF-system

The following page is derived from a project report written by Rafael Flock and Luca Sardi for the course Modelling and Simulation in Structural Mechanics in the winter term 2019/20.

FE-Discretization of the Beam-like Resonator

Reference Model

The first and most straightforward way to discretize the resonator is to approximate the resonator by several Euler-Bernoulli beam elements with different cross section properties for the different elements. Note that the cross section is constant along each element.

Discretized resonator

The thin slices of the inclined part and the cubical mass at the end of the beam have dimensions which wouldn't allow a modeling with Euler-Bernoulli beams at first glance. But since this end part of the resonator doesn't deform very much during the oscillations it is assumed that the Euler-Bernoulli theory is still valid.

Tip Mass Model

Alternatively the resonator can be modelled with a single Euler-Bernoulli element with an additional rigid mass at the free end. Through this approach the resonator can be discretized by only two DOFs.

The FE-element of the tip mass model

The FE mass matrix for the Bernoulli-beam element with DOFs number one and two fixed reads as follows:

(1) \begin{equation} \mathbf{M} = \frac{\mu l}{420} \, \begin{bmatrix} 156 &-22l \\ -22l &4l^2 \end{bmatrix} \end{equation}

According to the Principle of Hamilton the prismatic mass at the tip of the beam introduces additional kinetic energy to the system. With m being the mass of the prisma and I_\Theta being its rotational inertia the additional energy is

(2) \begin{equation} \Delta T = \frac{1}{2} \, \dot{w}_2^2 \, m + \frac{1}{2} \dot{\varphi}_2^2 \, I_\Theta \end{equation}

Hence the Lagrange function is

(3) \begin{equation} L = \Delta T - \Delta U = \Delta T \end{equation}

In order to get the differential equation of the two active DOFs, variational calculus can be applied on the Lagrange function. This is done by using Euler's differential equation for each DOF. The resulting equations read

(4) \begin{align} m \, w_2 + 0 \, \varphi_2 = 0 \\ 0 \, w_2 + I_\Theta \, \varphi_2 = 0 \end{align}

Since we are considering a linear system the principle of superposition is valid and we can add m and I_\Theta to the corresponding entries of the mass matrix (1) in order to obtain the modified mass matrix \textbf{M}_{tm} of the tip mass model.


(5) \begin{equation} \mathbf{M}_{tm} = \frac{\mu l_{tm}}{420} \, \begin{bmatrix} 156 &-22l_{tm} \\ -22l_{tm} &4l_{tm}^2 \end{bmatrix} + \begin{bmatrix} m &0 \\ 0 &I_\Theta \end{bmatrix} \end{equation}

The mass part represents a prismatic solid, therefore the formula for the rotational inertia is given as:

(6) \begin{equation} I_\Theta= m \, \left( \frac{I_{yy}} {A} + \frac{t_{m}^2} {12} \right) \,\, \text{with} \,\, m=\rho \, A \, t_{m} \end{equation}

A is the area of the tip mass surface and I_{yy} the moment of inertia around the y-axis w.r.t. the center of gravity

Topview on the tip massNote that l_{tm} is not equal to the actual length  of the resonator but represents the distance between the fixed end and the center of gravity of the mass of the tip mass model. Hence the actual length l_{1} of the resonator can be obtained as:

\begin{equation} l_1=l_{tm}-x_G \label{eq:act_l1} \end{equation}

Tuning of the different Models

Tuning of the reference and the tip mass model

The approach for the reference model and the tip mass model is to fix all properties of the resonator related to geometry and material except for l_{1}. It is also possible to choose a different free parameter, e.g. l_{3}l_{1} however proved to give the best match between the two models. With l_{1} representing a free variable the stiffness matrix \textbf{K} and mass matrix \textbf{M} can be written analytically as functions depending on l_{1}. The target frequency to which the resonator has to be tuned can be written as circular frequency \omega. Solving the relation

(7) \begin{equation} \det\left(\mathbf{K}(l_1) -\omega^2 \, \mathbf{M}(l_1)\right) = 0 \end{equation}

which normally is used to determine the eigen frequencies of a given system, different solutions for l_{1} can be found. Imaginary or also negative real solutions for l_{1} can be neglected since l_{1}, being the length of an actual beam has to be real and positive. Among the remaining solutions, the one that leads to the first eigen frequency of the resonator being equal to the target frequency is the correct one.

This procedure is implemented in MATLAB by defining equation (7) symbolically as function of l_{1} for the given resonator properties and target frequency \omega. The different solutions of equation (7) are found analytically with the command solve. To finally find the right solution among the different found solutions a modal analysis of the resonator for the different real and positive solutions of l_{1} is carried out.

Tuning of the SDOF Model

Neglecting damping, the two parameters by which a SDOF system is described are the mass m and the spring stiffness k. The mass can thereby be derived directly from the geometric properties of the resonator. An other option to obtain the mass m is to sum up the entries of the mass matrix \textbf{M} of the reference model which are related to the vertical nodal forces due to a vertical motion of the beam.

(8) \begin{equation} m = \sum_{i=0}^{\frac{n_{DOF}}{2}-1} \sum_{j=0}^{\frac{n_{DOF}}{2}-1} \mathbf{M}_{(2i+1,2j+1)} \label{eq:sum_M} \end{equation}

The two approaches give the same results for m if the model mass matrix \textbf{M} represents the exact geometry. But since the skewed part of the mass is discretized by piecewise straight Euler-Bernoulli beams there is a small difference between the two results. The second approach, which gives the exact mass of the reference model, but not of the originally given geometry, is applied.

Since the eigen frequency of the SDOF system has to be the same as the target frequency the relation

(9) \begin{equation} \omega^2 = \frac{k}{m} \end{equation}

is used to determine the appropriate spring stiffness k of the SDOF model.

Ideal hysteretic Damping

Damping is introduced by modeling Young's modulus as a complex quantity:

(10) \begin{equation} \hat{E} = (1+\eta \, i) \, E'' \end{equation}

The parameter \eta is called loss factor and describes the ratio between the loss and the storage (elastic) modulus which are denoted here by E' and E'', respectively:

(11) \begin{equation} \eta = \frac{E'}{E''} \end{equation}

It must not be confused with the percentage of critical damping D, also denoted as \xi . The relationship is eta=2D .

Since the Young's modulus appears as a factor in front of the stiffness matrix of Euler-Bernoulli beam elements the damping matrix and the stiffness matrix of the dynamic stiffness can be combined in a complex stiffness matrix:

(12) \begin{align} \mathbf{D} &= -\Omega^2 \, \mathbf{M} + i \, \Omega \, \mathbf{C} + \mathbf{K} \nonumber \\ &= -\Omega^2 \, \mathbf{M} + i \, \eta \, \mathbf{K} + \mathbf{K} \nonumber \\ &= -\Omega^2 \, \mathbf{M} + \mathbf{K} \, (\eta \, i + 1) \end{align}

Hence, the damping matrix is

(13) \begin{equation} \mathbf{C} = \frac{\eta}{\Omega} \, \mathbf{K} \end{equation}

This is a special case of Rayleigh-damping where \textbf{C}=\alpha\textbf{M}+\beta \textbf{K}, with \alpha=0 and \beta=\frac{\eta}{\Omega}. As a consequence the eigen vectors and eigen values derived in the preceding sections with \mathbf{K} and \mathbf{M} don't change and the tuning approach is still valid. This means that, in order to include the damping of the resonators into the resonator models of the unit cells, it is sufficient to adjust their generalized stiffness by the factor (1+i \eta). Since the damping model and \eta is valid for the host structure as well and the excitation is harmonic, not more than multiplying the assembled stiffness matrix of the complete system with that factor is necessary.

Discussion and Comparison of Results

Parameters

The parameters for the following results are given in table. Note that w_{m} is given implicitly by the maximal angle of fabrication (30 degree), w_{B} and l_{2}.

Parameters of the test configuration

WFEM Analysis of the tuned Models

For the following examples the target frequency is the first eigen frequency of the host structure, i.e. 110.63 \text{Hz}

Host structure for the validation of the different models

All models are tuned according to the methods described in the preceding section. The resulting specifications of the models are given in the following table. The relative mass is given w.r.t to the mass of one unit cell.

Results of the Tuning

A WFEM analysis was carried out for three unit cells, equipped with the reference, the tip mass and the SDOF model, respectively. The tolerance measures were 1.00\text{cm} for the decay and 0.001\text{rad} for the phase change. The corresponding plots of the decays per unit length are given in the figure below. Evanescent shares are marked in purple, propagating shares are marked in red and the black dots show the frequencies where the unit cells are damped. The frequency for which all three cells show damping behaviour is the target frequency, i.e. 110.63 \text{Hz}. In the plots for the reference and the tip mass model a second bandwidth with damping characteristics is visible: Those are the second eigenfrequencies (1973.3 \text{Hz} and 1311.6 \text{Hz}). They do obviously not match as well as the first eigen frequencies.

Decay per unit length

All three unit cells have almost the same stop bands at the target frequency. In order to show this, a more detailed WFEM analysis was carried out around the target frequency (with 0.001 \text{Hz} - steps) and the first and last frequency noted for which attenuating waves occurred.

Bounds of the stop bands at the target frequency

Host structure equipped with the Reference Model

The following figure shows the first two eigen modes of the host structure ('hs'), i.e. the pinned beam.

The first two eigen modes of the host structureBy adding ten equally distributed resonators (modelled with the reference model, 'rf') to the structure, each tuned to the first eigen frequency of the host structure (at 110.63 \text{Hz}), eleven new eigen modes appear: Nine of them have more or less the former first eigen frequency (only one shown in the figure) and two of them are located in the vicinity below (109.26 \text{Hz}) and above (112.02 \text{Hz}) the target frequency. The next eigen mode at 442.53 \text{Hz} is more or less the same as before, i.e. the former second eigen mode. All eigen modes are normalized w.r.t the mass matrix.

Eigen modes of the host structure equipped with the reference modelThe host structure is damped very well at its former first eigen frequency. The displacements of the beam-like resonators are very high in the given graphs but since we are investigating the eigen modes here, they are scalable and represent only the relation to the host structure which deforms much less. The new eigen modes below and above the target frequency show quite large displacements, consequently causing the frequency response in these regions to worsen as well.

Comparison of Responses of all models

The following figure shows the responses of the host structure equipped with the reference model, with the SDOF model ('SDOF') and with the tip mass model ('tm'), respectively, for selected frequencies. All resonator models show almost the same displacements. As already pointed out in the previous section one can see that the resonators have only a damping effect at the target frequency and its immediate vicinity (111 \text{Hz}). Below and above this frequency the displacement of the host structure gets even worse (109\text{Hz} and 112 \text{Hz}). At all other frequencies above and below this range almost no effect of the resonators is visible.


In the left figure one can see the maximal absolute displacement of the center of the host structure for all the three models w.r.t. the excitation frequency. The two new peaks of the damped structure are clearly visible. In between the peaks (109.26 \text{Hz} and 112.02 \text{Hz}) only within a small bandwidth good damping is achieved. Below the lower and above the upper peak the host structure behaves more or less the way it would without resonators.

The SDOF and the tip mass model behave equally well w.r.t. to their amplitudes and the corresponding frequency, whereas the tip mass model shows a slightly better but practically negligible fit to the reference model (see right figure and the table below).

Frequency responses

Zoom into left figure

As can be seen in the table below the maximal displacement of the structure is reduced from approximately 375 \text{mm} to 190\text{mm}. On top of that, the first eigen frequency is successfully damped. But the huge disadvantage is the appearance of the new peaks in the vicinity of the target frequency which immediately leads to high displacements even if the excitation frequency deviates only weakly from the target frequency.

Differences w.r.t amplitudes and corresponding frequencies

Results for different target frequencies and higher relative mass

To further validate the different modelling approaches the complete procedure of tuning, assembling and analysis of the frequency response is carried out for different lengths of the quadratic cross section of the host structure. The side length is varied from 0.01 \text{m} to 0.05 \text{m}. With this change the first eigen frequency changes and now varies from approximately 22 \text{Hz} to 111\text{Hz}. The larger the cross section becomes, the more the first eigen frequency increases. Furthermore all dimensions of the resonators are doubled in order to assure a higher relative mass in combination with the smaller cross sections of the host structure. The force amplitude is adapted to the smaller cross section of the host structure, meaning it is decreased to 10.0 \text{N}. The changed parameters are given for convenience in the following table:

Changed parameters

The figure below shows the maximal absolute displacement in the center of the beam and the corresponding relative reduction compared to the undamped host structure. The maximal displacements of the systems equipped with the resonator models correspond to the lower of the two due to the damping emerging peaks which were described in the previous section.

All models still behave equally well. The displacement in the first graph decreases with higher eigen frequencies mainly since the cross section becomes larger and thus the host structure stiffer. The effect of the resonators is visible in the second graph: A displacement reduction of approximately 40 to 48 \% is achieved for the given eigen frequencies.

Absolute displacement and displacement reduction for different target frequencies


In the figure in the first graph one can see the resulting lengths  l_{1} from the tuning for the different eigen frequencies. They decrease with increasing stiffness and hence increasing eigen frequency of the host structure. Analogously the relative mass decreases in the second graph since on the one hand side l_{1}gets smaller as already noted and on the other hand side the system's mass gets larger due to the larger cross sections.

Tuned length l1 and corresponding rel. mass

References

Kausel, Eduardo. Advanced structural dynamics. Cambridge University Press, 2017.

Wikipedia contributors. Dynamic modulus —Wikipedia, the free encyclopedia,

2019. [Online; accessed 27-February-2020].