Frequency filter for elastic bending waves: Poincaré map method and experiment

In this work, the Poincaré map numerical method was successfully developed to solve the fourth-order differential equation that describes the flexural vibrations of a beam, within the Timoshenko beam theory. The Euler-Bernoulli continuity conditions were considered, which are valid for frequencies smaller than the critical frequency. As an example, this method was used to design a complex elastic structure, characterized by a flexural frequency spectrum with a broad band gap. Such structure consists of two coupled phononic crystals, which were designed with filling factor values in such a way that in their bending frequency spectra, an allowed band of the first part, overlaps with a band gap of the second one and vice versa. The resulting composed system has a much wider effective gap than its original components, between 4 and 10.5 kHz. This system works as an elastic bending wave filter. Finally, these three structured elastic systems were constructed, and characterized by the acoustic resonance spectroscopy technique. The natural flexural frequencies as well as the corresponding wave amplitudes of each structured beam were measured. The experimental measurements show excellent agreement with the numerical simulation.


INTRODUCTION
The study of acoustic and mechanical periodic systems have matured into a very active research field since physicist established an analogy to solid state in the early 1990's [1]. Ever since, the so-called phononic crystals have become a setting for investigating novel phenomena related to topological phases [2]. Acoustic and elastic waves propagating through a periodic arrangement of structures show frequency regimes in which the waves stop propagating, this phenomena is equivalent to the electronic bandgaps in semiconductors [3].
On the other hand, some of the most common methods to make theoretical calculations of spectral band engineering in locally periodic structures are: the finite element method [30,31], spectral element method [32], the transfer matrix and modified transfer matrix methods [33]. One particular method to simulate vibrations in continuous media based on the optimization of beam systems are the spring-mass arrangement [34][35][36][37] However, it has been shown that the Poincare map is faster, more capable and stable than the transfer matrix method [38], particularly for compressional and torsional waves, and it has not been implemented for flexural waves to the best of our knowledge.
The Poincaré map method has also been used as an effective tool in damage assessment for the problem of vibration health monitoring [39], to verify aperiodic behavior of ultrasonic waves due to load-induced cracks in concrete materials [40], and it has been proposed as an accurate method in structural damage detection [41]. Nevertheless, it has never been used to investigate elastic bending waves in the design of phononic crystals, which are phenomena described by a fourth-order differential equation. For this reason, designing and manufacturing phononic crystals which work with this type of waves is inherently more difficult as this study shows.
In this work, the Poincaré map method, successfully developed to compute the bending spectrum in arbitrary one-dimensional structured systems, is reported. In the next section, this numerical method is developed and used to calculate the normal-mode frequencies of a one-dimensional (1D) finite phononic crystal with different filling fractions. In Section 3, the Poincaré map method is applied to calculate a structured system consisting of two coupled finite phononic crystals. This composite structure and two finite phononic crystals were constructed. Each of these systems were machined from a single aluminum circular beam. Numerical and experimental results are compared, and show a very good agreement.

A POINCARÉ MAP FOR BENDING WAVES
Some of the methods commonly used to calculate bending vibrations of structured beams are the finite element method and the transfer matrix method as mentioned above. The former, has the disadvantage that the results come combined with other kinds of vibrations, such as torsional and compressional ones, and it is complicated to separate the different type of vibrations from the one of interest [42]. The latter, can only be applied to structures with a very small number of pieces, since it becomes slow and unstable when there is a large number of parts [43].
The Poincaré map method was originally developed for the study of 1D heterostructures [44]. After that, it was applied to calculate the torsional and compressional spectra of vibrations in beams [29]. This extension was natural since both the quantum and the classical wave equations are second order equations and the main differences are the boundary conditions. A completely different situation occurs with bending vibrations since they are ruled by the Timoshenko beam theory, which is described by a fourth-order differential equation. Thus the extension of the Poincaré map method for bending waves is much more complex.
The theory which describes flexural vibrations of beams is known as the Timoshenko beam theory [45]. Apart from the transverse displacement ξ , which appears in the Euler-Bernoulli (E-B) beam theory, another degree of freedom exists, the angular coordinate [46]. The system of coupled equations, which satisfies ξ and , is given in the Appendix. This system of equations can be written as a single equivalent fourth-order equation for the transverse displacement, known as the Timoshenko beam equation. For a uniform beam of cross-section area S, shear modulus G, Young modulus E and density ρ, this equation is given by An identical equation holds for the angle . Here κ is the Timoshenko shear coefficient, and I is the second area moment with respect to the axis perpendicular to the rod axis. Even though it seems paradoxical that only the transverse displacement ξ appears in the previous equation, the angle of rotation can be written as a combination of the first and third spatial derivatives of ξ (see Eq. (59) in the Appendix). Imposing the normal mode condition ξ (z, t) = χ (z)e iωt , where i = √ −1 and ω = 2πf, being f the frequency, the Timoshenko equation can be reduced to where and is the critical frequency. For a uniform rod, the amplitude χ (z), solution of the last equation, reads as where being [j/2] the largest integer less or equal than j/2. Notice that the factors k −3 j multiplying the coefficients A j are arbitrary since they can be absorbed in the coefficient A j .
Let us consider a rod of length L composed by n coaxial cylinders of cross-sectional areas S i , as shown in Fig. 1. The amplitude χ i (z) in the i-th cylinder is where z i ≤ z ≤ z i+1 , i = 1,…, n and k j, i with j = 1,…, 4 is given by: where and is the critical frequency.
Here ω is the angular frequency, G i and E i are the shear and Young's modulus, respectively, ρ i is the mass density, S i is the cross-sectional area, κ i is the Timoshenko shear coefficient and I i is the second moment of area respect to the perpendicular axis of the rod axis. All the above expressions correspond to the i-th cylinder. Conditions of continuity at the points where the cylinders of different radii join each other are as follows: continuity of the displacement, slope, bending moment and shear force: The amplitude χ i (z) and its derivatives depend on the physical and geometrical properties of i-th cylinder through Eqs. (9)- (12). It should be mentioned that these continuity conditions are approximate since the beam is three-dimensional (3D) whereas the model is 1D [29]. Also, only the regime of low frequencies is explored in the experiment, compared to the lower critical frequency given by Eq. (12). Thus, the E-B continuity conditions are used since in the limit ω ω c,i ,i = 1,…, n, (see Appendix), the more accurate Timoshenko continuity conditions are reduced to the E-B ones [47]. Also, when the second equation in Eq. (13) vanishes, the continuity conditions yield a quadratic function, which is a good approximation to the exact wave amplitudes. This effect can be seen in the extensive studies of bending wave amplitudes in structured beams [48]. Let's consider the vector where the superscript T means transpose. Using Eq. (8), the last equation can be written as where and is a 4 × 4 diagonal matrix. By applying the boundary conditions, Eq. (13) into Eq. (15), it is obtained where ζ i ≡ ζ i (z i ). K i is a 4 × 4 matrix given by Rewriting the matrix K i as composed by a set of 2 × 2 matrices thus, by defining two vectors and Eq. (18) can be rewritten as two separate equations as follows and Eq. (23) is solved for η i , Shifting the index i → i + 1 in this equation it is obtained Substituting Eqs. (25) and (26) into Eq. (24) one gets an equation exclusively in terms of two-component vectors φ By inspection of Eq. (9), it can be seen that k 1, i = −k 3, i and k 2,i = −k 4,i . By using this, simplified equations can be written and where Additionally, these matrices satisfy the following identity which allows to rewrite Eq. (27) in a more compact form then, the trivial result is Applying the boundary condition for the free end at z 1 = 0, one gets χ 1 (z 1 ) = χ 1 (z 1 ) = 0. Eq. (22) yields η 1 = 0 and therefore from Eq. (23) it is obtained and hence L 2 = O 1 . Substituting L 1 and L 2 into Eq. (36) with i = 1 the result is thus, Eq. (36) allows to know L i+2 in a recursive way The normal-mode frequencies are obtained when the free-end boundary condition, χ n+1 = χ n+1 = 0, is satisfied. Substituting this into Eq. (23), with i = n, yields or equivalently Notice that this is a condition for a 2 × 2 matrix equation, unlike the transfer matrix approach, which deals with 4 × 4 matrices. Although the Poincaré map approach needs more arithmetical operations, it is numerically more stable than the transfer matrix method and can be used for a larger number of elements in the structured beam.

MECHANICAL FILTER FOR BENDING WAVES
Two locally periodic beams were studied to illustrate the advantage of the Poincaré map method developed in the previous section. In order to find the optimal parameters for both beams, Eq.   Beams A and B were chosen in such a way that the gap in the spectrum of beam A overlaps with the passband of the spectrum of beam B and vice versa. Beam C is built as the composition of beams A and B and presents a larger bandgap since the allowed modes of a band of the first phononic crystal will be blocked by the gap of the second phononic crystal and vice versa, as shown in Fig. 2(c). This is, beam C works as a wide frequency filter for bending waves. A similar 1D periodic aluminium beam for studying the change in topological phases has been recently published [21]. The flexural wave spectrum for the beam of Fig. 2, as a function of the filling fraction F = / , is given in Fig. 3. Here is the width of the notch and the size of the unit cell. It can be seen at low filling fractions that two bands and one gap appear. The lowest band starts at frequency zero due to the free-free boundary conditions and the gap becomes smaller as the value of the filling fraction increases and almost disappears around F = 0.6.
At low filling fractions, within the interval of the plot, two bands and one gap appear. Due to the free-free boundary conditions, the lowest band starts when the value of the frequency is zero. Then, the gap becomes smaller as the value of the filling fraction increases and almost closes at F ≈ 0.6. At higher frequencies, there is a second allowed band. The case F = 0 corresponds to = 0 and then the beam is uniform with diameter d . Conversely, F = 1 corresponds to = and the beam is again uniform with a diameter d . It is possible to identify two optimal values of the filling fraction, F A = a / = 0.144 and F B = b / = 0.499, which characterize the two first periodic designs (systems A and B in Fig. 2). These values were selected because the first gap ends in a different frequency for each beam, and the second band starts at different frequencies too.
Systems A and B have their first gap within 5 f 12 kHz and 6 f 7 kHz, respectively. As shown in Fig. 2, system C is a combination of two locally periodic structures of beams A and B. It exhibits a wider gap in the interval 5 f 12 kHz since the composite structure does not transmit frequency signals in the gap of both structures. Figure 2 shows drawings of beams A, B and C, presented as (a), (b) and (c), respectively. These beams were machined, and then characterized by using the non-contact acoustic resonance spectroscopy technique [49]. A diagram of the experimental setup used to characterize the spectra of the structured beams is shown in Fig. 4. It consists of a vector network analyzer (VNA) Anritsu MS4630B (1); a high fidelity audio amplifier (Cerwin-Vega CV900) (2); two electromagnetic-acoustic transducers (EMATs) (3) and (5); a workstation (6); and the aluminum beam under study (4). The VNA generates a harmonic signal of frequency f, which is sent to the audio amplifier. The output of the amplifier, in turn, is sent to a second EMAT exciter. The bending waves, at one end of the rod, are induced without contact by the EMAT exciter (3). Simultaneously, at the other end of the rod, the EMAT (5) measures the vibratory response of the beam to the excitation. The detected signal is then registered by the VNA to be compared with the monochromatic original signal. All data are sent to a computer for analysis. Finally, the frequency of the exciter signal is increased by f and the cycle starts again; then a given frequency interval is scanned.

ACOUSTIC RESONANCE SPECTROSCOPY RESULTS AND COMPARISON WITH THEORY
The bending waves are excited selectively on the rods using a particular configuration of the EMATs. The electromagneticacoustic transducers consist essentially of one coil and magnet, which are placed in particular configurations depending on the vibration type. Since the transducers operate through eddy currents, they have to be located in the vicinity of a paramagnetic metal, this is the reason why the beams used here are made of aluminum. A detailed explanation about the operation principles of these transducers can be found elsewhere [49,50]. Both EMATs configuration used to excite and to detect bending vibrations in the structured beams is illustrated in the inset of Fig. 4. The coil  Fig. 2. The left (cyan) column shows the results from the Poincare map method, the middle (blue) column shows the results of the finite element method-COMSOL Multiphysics-and the right (purple) column shows the data from experiment. The levels highlighted in pink of system C do not contribute to the transport; they also appear in the Poincare map method results. The horizontal dashed lines indicate the experimental wide band gap of C, which is very similar to the system A. In the right part, the first band is shown and the begging of the second one can be seen. and magnet are coaxial and their magnetic dipole is perpendicular to the rod axis.
Figs 5-7 show the measured bending spectra for beams A, B and C in different frequency intervals, presented as (a), (b) and (c), respectively. In Fig. 5, the spectra of the three beams from 400 Hz to 4 kHz are shown. The sharp peaks correspond to the normal-mode frequencies of each structured beam. The plots at the right side show zooms of the end of the spectra. Since the beams have cylindrical cross-section, the bending vibrations in these beams should be degenerated. However, because of manufacturing defects, the degeneration is broken and pairs of levels can be observed. These pairs of levels-which will be referred as doublets for practical purposes-are the experimental fin-gerprint of the bending vibrations in cylindrical (and squared) beams and can be used to distinguish flexural vibrations from torsional and compressional ones [49]. A small offset can also be observed at the right panels and in the subsequent figures; this is due to the impedance between the coils of both EMATs; this gives rise to Fano resonances [51].
The measured spectra at higher frequency interval of beams A, B and C are shown in Fig. 6 in the interval [4750-8750 Hz]. At the left part, zooms of the end of the spectra are shown. It can be seen that some unwanted resonances appear in beams A and C. These are associated to torsional or compressional vibrations. Besides, the spectrum of beam B shows the doublets of the second band that starts at 4975 Hz. Based on Fig. 3, a gap is expected for both beams A and C in this frequency range. In the case of beam B, the end of the first band and the beginning of the second one can be observed. The typical doublets of bending vibrations can be noticed in the spectrum of beam B, then the second band starts at ≈6 kHz. The spectrum shows different intensity for each resonance of the doublet since the intensity depends on the orientations of the EMATs respect to the principal axis of the rod. Figure 7 shows the experimental spectra of the beams in the range of frequencies from 10 to 15 kHz. It can be noticed that beams A and C show the beginning of their second band at ≈10.5 kHz. Once again, this is consistent with the numerical results of Fig. 3 for the filling fraction of beam B. Finally, beam C works as a filter of bending waves, since there exist a band if beams A and B have a band in a frequency interval. A gap appears in beam C in the frequency intervals in which either beam A or B have a band gap. The phenomenology of Figs 5-7 is illustrated in Fig. 8, in which a comparison between numerical and experimental bending spectra resulting of the three structured beams is shown. The finite element method results are also given for comparison. The spectra on the left correspond to the structured beam A of Fig. 2. In the middle, the spectra of beam B is shown and the last one at the right are those of the com-  ) and (c) are localized at the right part of the structured beam, therefore do not contribute to wave transport. The associated frequency levels to these modes lie inside the measured wide band gap and appear in the FEM and Poincare map method spectra (pink levels in Fig. 8).
posed beam C. The designed wide gap of the composite structured rod (C) is almost the same width as the one in beam A. This is because of the corresponding gap of beam B is small and is completely contained in the gap of beam A, in this frequency range. Besides, at higher frequencies above 20 kHz, a wider gap is also expected in beam C, as it can be seen in Fig. 3. Some normal modes computed by using finite elements are given in Fig. 9; modes (a) and (d) are extended whereas modes (b) and (c) are localized.

CONCLUSIONS
In this work, the Poincaré map method was successfully developed to calculate elastic bending vibrations in structured beams. It is the first time that this method has been applied to a fourthorder differential equation that describes bending waves. The Poincaré map method allows the calculation of the bending frequency spectra in arbitrary structured long beams, for which the conventional transfer matrix method does not work due to instability conditions. The capability of the numerical method presented here was tested by engineering a wide bandgap mechanical filter composed by two phononic crystals. The giant bending-waves bandgap of the engineered system ranges from 4 to 10.5 kHz. The width of this gap is the result of the overlap of the corresponding gaps of each phononic crystal. The structured beam was fabricated in aluminum and its flexural frequency spectrum was measured by using acoustic resonant spectroscopy. The experimental results show an excellent agreement with the theoretical predictions. Giant-gap mechanical filters have potential applications as vibration reducers in buildings, walls, bridges, machines and vehicles in general. Additionally, scaled versions could be used as selective mechanical sensors or transducers.

APPENDIX : CONTINUIT Y CONDITIONS FOR THE TIMOSHENKO BE A M THEORY
In this section, the Timoshenko continuity conditions are derived by following a procedure similar to the one given in [52] to obtain the boundary conditions for a free uniform beam. Then, in the limit of low frequencies compared with the critical frequency, the E-B continuity conditions in Eq. (13) are obtained. The boundary conditions for a change of cylinder ends are given by the continuity of transverse displacement ξ and rotation angle and the continuity of bending moment and shear force, respectively. Here the subscript L(R) labels the cylinders at the left (right) of z i . For the time-independent variables ψ, and χ , where (z, t) = ψ(z)e iωt and ξ = χ (z)e iωt , the first two continuity conditions, Eqs (45) and (46), trivially become while Eqs (47) and (48) yield where the first one, Eq. (49), corresponds to the first continuity condition in Eq. (13), and no approximation is involved. Since the continuity conditions in Eqs (49)-(52) mix variables ψ and χ , they have to be reformulated only in terms of one or the other. In particular, we have chosen χ as the variable of our study both theoretically and experimentally. So the continuity condition in Eq. (50) involving only the angle ψ has to be reformulated in terms of χ . Let us consider the coupled equations (See Eqs (1) and (2) of reference [52]) equivalent to the Timoshenko equation For the time independent variables ψ and χ , one has Solving for dψ dz from Eq. (55) one obtains and deriving Eq. (57) one gets Substituting Eq. (58) into Eq. (56) and solving for ψ yields The continuity of the angle of rotation ψ L z=z i = ψ R z=z i then can be written in terms of the derivatives of χ as follows, When ω ω c , the angle of rotation ψ seems to be given by a linear combination of the first and third derivatives of the transverse displacement χ However, by using Eq. (6), the first derivative of χ is given by while the third derivative is In the limit ω ω c , the square of the wave number k 2 j , given by Eq. (7), is therefore, the first term in Eq. (61) is EI κGS which turns out negligible compared with the second term of Eq. (61), dχ dz , and given by Eq. (62). In Eq. (65), we have replaced E G by 2(1 + ν), where ν is the Poisson's ratio. Therefore, in the limit ω ω c , Eq. (61) reduces to the Euler-Bernoulli continuity condition, and Eq. (50) reformulated in terms of χ , reads as corresponding to the second continuity condition in Eq. (13). Notice that this equation is no valid when dχ dz | z=z i = 0. In this case, the wave amplitude yields a quadratic function with curvature given by the second derivative of the transverse displacement, i.e. by the bending moment. This effect can be observed in the extensive wave amplitude measurements of [48], in which the discontinuty of the slope in the Timoshenko beam theory formulation appears as a quadratic function in the experiment. A similar behavior can be observed in the older experiments of [49].
By substituting Eq. (57) into the continuity of bending moment condition Eq. (51), one gets where ω c,P = κ P G P S P ρ P I P is the critical frequency for beam P = L, R. For low frequencies, compared with the critical frequency, one has ω ω c, P , and then the third E-B continuity condition in Eq. (13) for the second derivative is obtained Note that if a similar argument as the one used in the derivation of the continuity condition in Eq. (66) is applied to the first terms in the square brackets of Eq. (67), it will introduce a factor ω/ω c, P . Nevertheless, the dominant terms would still be the first terms in the square brackets in the limit ω ω c, P . Solving for the term in parenthesis in Eq. (56) and substituting into the continuity of shear force condition in Eq. (52), one obtains Once again, for frequencies much smaller than the critical frequency, ω ω c, P , the E-B continuity condition for the third derivative is obtained corresponding to the fourth continuity condition in Eq. (13).