Numerical Analysis of Vibration Attenuation and Bandgaps in Radially Periodic Plates

Periodic configuration of mechanical and civil structures has shown great potential for noise and vibration reduction. However, the use of Cartesian coordinates in studying periodicity effects in elastic structures overlooks the benefits of radially periodic configurations when dealing with wave propagation in large flexible plates disturbed by a small source area. This paper presents an easy-to-use numerical approach to predicting bandgap characteristics in polar coordinates. To demonstrate the vibration-attenuation effect, we consider a circular radially periodic plate model. We use an adapted Wave Finite-Element method in numerical experiments to demonstrate the existence of the attenuation effect. To verify the numerical results, we apply an adapted Floquet theory to polar coordinates. Our findings indicate that theoretical and numerical results are in excellent agreement considering a new parameter that introduces the distance from the origin. The adapted Wave Finite-Element approach and Floquet theory presented here demonstrate their potential to model more complex structures in polar coordinates.


Introduction
Periodic structures have the properties to strongly attenuate vibration transmission in specific frequency ranges, called stop-bands [1].This periodicity effect is often called Bragg scattering effect with reference to the propagation of X-rays in crystalline materials investigated in 1912 by William Lawrence Bragg and William Henry Bragg.Due to its potential applications in noise and vibration harshness mitigation in various engineering areas, modelling periodic structures have become an active and prolific research subject in structural dynamics.The interest has grown with the new technological possibility of creating multi-scale periodic structures with local resonance phenomena, which has made the design of these structures very attractive for new metamaterials [2][3][4].
Floquet's theorem, applied in 1946 by Brillouin to solve the wave equation [1], is a powerful approach to studying waveguide properties of periodic structures from the knowledge of the dynamic properties of a single period (cell).Since this theorem heavily relies on the translational invariance of the problem formulation, periodicity effects have been studied by its means mainly in Cartesian coordinates, while polar periodicity has very seldom been examined, for example [5][6][7].However, many engineering applications deal with structures whose dynamics can be conveniently described in polar coordinates, e.g.flexible membranes or plates of a circular shape, built-up structures with beam/ plate connections, large flexible structures locally excited by alternative machines, and so on.In all these cases, radial periodicity can be utilised to reduce vibration transmission, and the periodic inserts can function as a vibroacoustic filter when the dominant excitation frequency falls within a stopband.The use of radial periodicity for vibration reduction was, for example, investigated in [8] where it was shown that a periodic structure, arranged in a stepwise shape along the radial direction, was able to attenuate vibration due to its low-frequency bandgap characteristics.Radial periodic structures have also been shown to be effective in protecting constructions from seismic waves; as an example, in [9], a new type of seismic metamaterial with radial periodicity was studied, demonstrating its effectiveness, compared to traditional seismic metamaterials, in attenuating seismic Lamb and surface waves.
Although a mathematically rigorous reformulation of Floquet's theorem for polar coordinates is problematic, its far-field analogue, known as the theory of Bragg fibre, has been formulated in optics [10][11][12].In [7], the Floquet propagator has been introduced as a function of the radial coordinate, and several heuristically chosen approximations of this dependence have been compared with each other and verified via direct numerical integration.The recent paper [13] has addressed the formulation of a propagator, which is equally applicable to Cartesian and polar coordinates.Then its independence upon the axial Cartesian coordinate naturally emerges, and its dependence upon the radial polar coordinate becomes uniquely defined in agreement with findings presented in [7].
The wave and finite-element (WFE) method has now become a classical tool to analyse properties of periodic and authentically uniform waveguides of arbitrary complexity in Cartesian coordinates.It has been presented in many papers and studies, and several benchmark cases have shown its practical applicability to periodic and continuous structures, e.g.[14].The method models a single period, or cell, using standard finite-element (FE) analysis, and it applies Floquet/ Bloch theory [1] to evaluate the behaviour of the whole structure in terms of wave propagation.Recently, the WFE has been generalised to predict free and forced wave propagation in radially periodic plates [15,16].In particular, in [15], the response to harmonic excitation was studied close to the excitation zone and in the low-frequency range (in the first pass-band), while in [16], the prediction of vibration attenuation was verified by comparison with standard finiteelement analysis of a finite plate.It is of worth mentioning here that stop-bands rigorously occur in infinite periodic structures, while attenuation bands are found in finite structures with the insertion of periodic cells.When the number of inserted periodic cells is large, it has been shown that the attenuation-and stop-bands are the same [17][18][19].
For consistency, the procedures outlined in [13] and [16] are presented to predict and verify stop-bands and vibration attenuation bands in radially periodic plates.The main advance from [13] and [16] is the generalisation of these procedures for an arbitrary circumferential wavenumber.The paper is organised as follows.In "Floquet Propagator in Polar Coordinates", the formulation of the Floquet propagator, equally applicable to Cartesian and Polar Coordinates, is presented, showing that it can be used to assess the WFE applicability for the identification of pass-bands in periodic plates in the same way as it has been done for membranes in [13].In "The Wave/Finite Element and Wave Amplitude Decay Prediction", the WFE approximation to model free and forced wave propagation for an arbitrary circumferential wavenumber is described.Only a few segments of a sector of the structure are discretised using conventional FE.Periodicity conditions are applied to develop a WFE polynomial eigenproblem whose solutions yield the wavenumber-frequency relation and wavemodes.These are further used to obtain information on the wave energy flowing and describe the structure's response to localised excitation sources."Numerical Examples" contains examples, which demonstrate excellent agreement between results obtained using the procedures outlined in [13] and [16] and their verification.

Floquet Propagator in Polar Coordinates
A schematic representation of a polar periodic plate is shown in Fig. 1.The figure shows a plate with repetition both in the circumferential direction, periodic angle , and in the radial direction, periodic length L. To avoid numerical errors and singularities, the first cell is shifted at an arbitrarily small distance R 0 from the origin of the polar coordinates.The internal radius R 0 is an important novel geometrical param- eter (in comparison with the problem formulation in Cartesian coordinates) of the problem.
In this paper, we consider only radially periodic plates.However, since vibrations of these plates feature a circumferential wavenumber m = 0, 1, 2, ... as cos(m ) , the set of equally spaced rays in Fig. 1 represents radial nodal lines of a vibration mode.In the framework of WFE modelling (see "Floquet Propagator in Polar Coordinates" for details), only one 'slice' of angle needs to be considered.Then, this angle may not obey the relation = m with m being an integer number.
As is well known for periodic structures in Cartesian coordinates, free waves-although purely time-harmonicbecome spatially composite, which is clearly seen in their Fourier decomposition [20,21].Conversion of the problem formulation and solution to polar coordinates implies the replacement of sine and cosine with Hankel functions.Then, a composite wave in the radially periodic membrane also has a discrete spectrum in Hankel decomposition rather than a single component.As substantiated in [7], formulation of Floquet propagator for such a membrane requires adjustment of its conventional Cartesian formulation using a "cylindrical wave correction factor", which is applied to the whole composite wave whatever its decomposition content.
To compare the stop-bands obtained with Floquet theory and the WFE method, we first briefly reproduce all required steps to determine the "cylindrical wave correction", or shape function, for a periodic plate.It is assumed that the dynamics of the plate is governed by the classical Kirchhoff plate equation shown in Eq. (1): where D is the bending stiffness expressed as with Young modulus E and Poisson ratio .In what follows, we use a time-harmonic state exp(−i t) and dimensionless frequency parameter Ω = L 2 √ h 1 ∕ D 1 and we assume that the structure has two parts with parameters

Cartesian Coordinates
For consistency, we reproduce here the classical solution of the one-dimensional "cylindrical bending" problem in Cartesian coordinates.General solutions of Eq. (1) for alternating parts have the form of Eq. ( 2): (1) For simplicity, we use dimensionless parameter in Eq. ( 2).To impose continuity of the structure, interfacial conditions are used as in Eq. ( 3): where M = w��(x) and Q = w ��� (x) are the bending moment and shear force, respectively.Using classical guidelines, we use periodicity conditions from Eq. ( 4): Obviously, the last two equations formulate periodicity of bending moments and shear forces.In Eq. ( 4), point x 0 (3) could be chosen arbitrarily due to the translational symmetry of operator in Eq. (1) in Cartesian coordinates.Furthermore, in Eq. ( 4), we introduce the Floquet propagator Λ , which is linked to the classical Bloch's propagation constant as Λ = exp (iK B ) .To find Λ , we substitute Eq. ( 2) into Eqs.( 3) and ( 4) and obtain a system of eight linear algebraic equations with respect to modal amplitudes.The resulting determinant is the fourth-order polynomial in Λ and has the form In the following subsection, we show briefly how the classical Cartesian formulation is modified to obtain the Floquet propagator in polar coordinates.

Polar Coordinates
In polar coordinates, the general solution of Eq. ( 1) acquires a form of a linear combination of Bessel functions of the same integer order, which defines the circumferential wavenumber.Therefore, solutions for alternating parts are In Eq. ( 6), H (j) m , j = 1, 2 are the Hankel functions of first and second kind of order m and I m , K m are modified Bessel functions of order m .Following the Cartesian guideline proposed in "Cartesian Coordinates", we use interfacial conditions to preserve the continuity of the system.However, the operator in Eq. (1) does not have translational symmetry in polar coordinates in contrast with the Cartesian case.Moreover, in contrast to the canonical formulation in Cartesian coordinates, periodicity conditions and propagator become dependent upon an initial hole size R 0 [7].Therefore, the classical mathematically rigorous formulation of the Floquet propagator Λ = exp (iK B ) with Bloch's constant becomes inapplicable and requires a modification.
Starting from interfacial conditions, all equations should contain R 0 as a parameter.Thus, the resulting inter- facial conditions have a different form from Eq. ( 3) and they are defined as shown in Eq. ( 7): (5) where and are the bending moment and shearing force, respectively.
Periodicity conditions and propagator should be functions of the initial hole size R 0 , and the modified periodic- ity conditions have the form The exact equation to find the propagator is obtained by substituting Eq. ( 6) in Eqs. ( 7) and (8).As substantiated in [7] and [13], the resulting determinant is a fourth-order polynomial in Λ: The difference from the Cartesian case is that coefficient a 0 (R 0 , Ω) is not constant.Therefore, we modify Floquet theory as shown in [13] with the assumption Λ(R 0 ) = S(R 0 )Λ , where S(R 0 ) is the shape function.The latter is obtained from equation a 0 ≡ 1 and has the form which matches the "cylindrical wave-correction" function for a homogeneous plate [13].Moreover, the same shape function appears for the homogeneous and periodic membrane.Then the propagator for a periodic plate in polar coordinates is defines the location of pass- bands.However, for a homogeneous membrane, the opposite case, , does not necessarily specify a stop-band zone if these propagators are computed in the immediate vicinity of the origin of coordinate system.Nevertheless, introduction of the shape function S R 0 permits to unify the formulation of propagators in Cartesian and polar coordinates and recovers the conventional criterion of stop-band formation "almost everywhere" in the Ω, R 0 plane.This result, known for a membrane model [13], is now generalised for the model of a plate, which, on top of propagating cylindrical waves, captures evanescent ones.

The Wave/Finite Element and Wave Amplitude Decay Prediction
Compared to the corresponding Cartesian periodic structure, in polar periodic structures, the wave amplitude attenuation is dependent upon the distance from the centre.Following the procedure in [16], this section presents how to model this decay for free and forced wave propagation using the WFE formulation.Continuous structures can be studied using this approach but assuming arbitrary periodic spatial discretisation [15].With reference to Fig. 1, if the plate is homogenous in the circumferential direction-radially periodic-the angle can be assumed to be arbitrarily small; if the plate is homogenous, non-periodic, both arbitrary small length L and angle are assumed.
A periodic slice of the plates is taken, and it is approximated by a finite number of piecewise rectangular waveguide as shown in Fig. 2. We assume that the periodic angle is small, so that the properties of the slice do not vary rapidly in the radial direction.The segments in Fig. 2, the unit cells of characteristics length L in the radial direction, are univocally identified by an integer number j = 1…n, representing the cell's radial position.
To apply the method, a unit cell j is taken and meshed using conventional FE analysis as in Fig. 2. The cell can be of complicated construction and with internal nodes.These are condensed and reduced to obtain a super-element, whose nodal degrees are ordered as Here R = T 3 T 4 T and L are the Degrees of Freedom (DOFs) at the left and right sides of the unit cell.The figure also shows that local coordinates can be rotated, by the angle , to model the desired curvature as in [22,23].The WFE method has been introduced in many papers, e.g.[14,15,[22][23][24][25], and it is not detailed here.For completeness, an overview of the method is given in the Appendix A. Compared to [16], the formulation in the Appendix is given for 1-and 2-dimensional axisymmetric structures.The latter allows the evaluation of circular plate modes separately.
Since the full libraries of commercial FE software can be used, the approach is not limited to thin circular or annular circular plates, but it can be equally applied to laminated, composite and thick plates using either 2D or 3D solid finite elements' formulation without more effort.Moreover, the inclusion of curvature, stress, temperature effects and damping in the model can be accommodated, if necessary, as described in [23].
For each segment j, complex dispersion curves, ( j , ) , nodal displacements j (wave modes) and nodal forces j , which occur under the passage of a wave, are obtained as described in Appendix A. Using these wavemodes as basis functions, the total nodal displacements and forces of a segment j can be expanded by the sum of a finite number of the positive and negative wavemodes so that where a is the wave amplitude and the + and -superscripts refer to positive and negative propagating waves [14,24].
An external disturbance, force or displacement, is applied to the left DOFs of the segment j = 1.Therefore, only positive going waves of amplitude + , propagating away from the excitation point, will be generated.As an example, if an external force + 1 + 1 = e is applied, the excited wave amplitudes will be obtained from (11) 9) in [14].To recover + 1 , instead of using pseudoinversion, it is advantageous to exploit the left eigenvectors of the WFE eigenvalue problem to avoid numerical issues.For the segment j, the left eigenvectors are denoted as Left eigenvectors are orthogonal to the wavemodes and can be normalised, = .Therefore, the excited wave amplitude for the first segment j = 1 is obtained by Once the forced wave amplitudes in the first segment are found, the wave amplitudes in the next segment can be predicted as follows.Since we are considering a lossless waveguide with slowly varying geometrical properties, we can assume that the power associated with each propagating wave is preserved along the waveguides.
In the corresponding Cartesian waveguide, the nodal displacements and forces between two adjacent cells are related by the periodicity conditions where + j (L) = diag exp(−i + j L) is the diagonal matrix with Bloch's propagation constant Using Eqs. ( 11)-( 15), the left nodal displacements and forces of the unit cell j + 1 are recovered from Using the left eigenvectors as in Eq. ( 14), the first estimation of the wave amplitudes in the j + 1 cell is obtained from In a uniform and lossless waveguide, the amplitudes of waves propagating in the structure do not change.This is not true in non-uniform waveguides, where the wave amplitudes are a function of the position.In particular, in our case, we are considering a lossless waveguide with slowly varying geometrical properties-due to the assumption of small angles -and we can assume that the power associated with each propagating wave is preserved along the waveguides.The time-average energy flow density in each segment j can be evaluated from where the superscript H denotes the Hermitian transpose.This energy flow is complex-valued [26]: its real part, Re(Π j ) , (13) being the average power; its imaginary part, Im(Π j ) , being the peak of the reactive power, which can be seen as the stored averaged energy density in analogy with circuit theory.The wave amplitudes change, due to the changes in the geometry, can be accommodated by requiring that the time-average energy flow density through the interfaces of each Cartesian segment is the same [25].Therefore, we deduce that the ratio between Re(Π j ) and Re(Π j+1 ) gives the wave amplitude decay (j+1) : Once the wave amplitude change is evaluated, the nodal displacements and forces of the next segment can be predicted as The same passages are repeated to obtain the displacements and nodal forces at the target radial distance.
In the stop-bands, where only evanescent and highly decaying waves exist, there is no significant energy flow but mainly reactive power and Notice that the same considerations hold if the waves are travelling towards the centre of the plate, but the wave amplitude growth will be given by the reciprocal of Eq. (18).

Numerical Examples
In this section, a numerical example of a plate with periodically varying thickness is presented.Figure 3 depicts a schematic representation of the plate with annuli of different ( 18) . The material properties are: Young's modulus 210 GPa, density 7850 kg m 3 and Poisson's ratio 0.3.The non-dimen- sional frequency Ω = L 2 √ h 1 ∕ D 1 is introduced, where , h 1 and D 1 are the density, thickness and bending stiffness of the corresponding homogeneous plate of thickness h 1 , and is the frequency in radiant.The geometrical characteristics of the unit cell are indicated in Fig. 3.
The nodal degrees of freedom of the cell in Fig. 3 are ordered as = T L T R T I T , where L, R and I are associated with the right, left and interior nodal degrees of freedom of the cell, with a similar expression for the nodal forces.In this case, the discretised equation of motion, = is partitioned as The internal DoFs of the unit cell must be reduced as described in [27].Dynamic condensation can be applied, viz.I = − −1 II IL L + IR R , or other more efficient con- densing approaches can be considered, e.g.[27].

Bandgap Analysis and Verification
Bandgaps can be clearly represented using dispersion diagrams.In particular, stop-bands occur when there are no propagating waves and, for lossless waveguides, the absolute value of the propagation is not equal to one, viz.
shape from Eq. ( 10), we can determine theoretical stopbands for different values of initial hole radius R 0 as shown in Fig. 4a.The Cartesian periodic Kirchhoff plate stop-bands are naturally recovered at the limit R 0 → ∞ .However, as seen in this figure, their convergence to the Cartesian stop-band is very fast-already at R 0 = L there is virtually no difference between locations of "polar" and "Cartesian" stop-band.For R 0 < L polar, Floquet propagator predicts location of the left border of this stop-band with reasonable accuracy, whereas the frequency of its right border is much higher than its Cartesian counterpart.It should be noticed that, as discussed in "Floquet Propagator in Polar Coordinates", conditions plate stop-bands when R 0 << L. To verify the positions of stop-bands predicted using Floquet propagator with cylindrical wave correction, a forcing problem should be solved; see [7].We do not address it here because we aim to compare two alternative approaches to the identification of stopbands.Figure 4b shows the first bandgap of the radial periodic plate obtained using the WFE approach and the Floquet theory described in "Floquet Propagator in Polar Coordinates" when the distance from the centre of the plate is R 0 = L.For values of R 0 > L, no significant differences are observed, while for values of R 0 < L the WFE model cannot provide accurate results due to the ratio between the sides of the unit cell, which result in FE mesh distortion errors in this example.
In radially periodic plates, cylindrical waves feature dependence upon angular coordinate theta as cos(m ) with a circumferential wavenumber m = 0, 1, 2, ... Figure 5a shows how the theoretical first stop-band varies for different circumferential modes and different R 0 .As the circumferential mode order increases, while R 0 remains small, the position of the first stop-band strongly shifts towards higher frequencies.However, all stop-bands converge to the Cartesian limit as the initial hole radius grows.It is hardly surprising because the cylindrical wave correction factor, Eq. ( 10), does not depend upon the circumferential wavenumber.Similarly to the case illustrated in Fig. 4, a forcing problem should be solved to validate these results.However, we are interested in the comparison of Floquet and WFE predictions.Figure 5b depicts the comparison of the first stop-bands for R 0 = L and different values of m using the Floquet theory in "Floquet Propagator in Polar Coordinates" and the WFE approach in "The Wave/ Finite Element and Wave Amplitude Decay Prediction".It can be seen that the differences in the results increase with the circumferential number.This is due to the differences in the models: Floquet theory in "Floquet Propagator in Polar Coordinates" assumes a classical theoretical approach for separating circumferential modes, while the WFE approximates the phase of change of a wave as it propagates around the circumference using Eq.(A.6) in the Appendix, see also [22].
Of course, plate periodicity generates many stop-bands.However, this subsection has been concerned only with the study of the location of the first stop-band in a radially periodic plate.The reason for this limitation is that, in the first place, we are interested in comparison of results, obtained using the alternative methods.On top of that, as frequency grows, wavelengths in constituents of a radially periodic plate become shorter and, therefore, "far-field" (where the difference between plane and cylindrical waves becomes invisible) emerges (i) at lower frequencies and (ii) closer to the origin of coordinates.This suggests that the location of high-order stop-bands in a radially periodic plate is very close to their Cartesian positions.Finally, advancing to high frequencies makes the applicability of the elementary Kirchhoff plate theory questionable, whereas a use of Timoshenko-Mindlin plate theory lies out of the scope of this paper.

Vibration Response and Attenuation to Local Harmonic Excitation
Vibration response of the radially periodic plate excited by a point force was analysed in [16] using the WFE approach, and the results were verified by comparison with those obtained by standard FEA of a correspondent finite periodic plate modelled in COMSOL Multiphysics ® .The FE model was a homogeneous circular plate of thickness h 1 and total radial dimension R = 10L , to which periodic annuli were inserted from R 0 = L.The plate was assumed to be free.We reproduce these results in Fig. 6, where the absolute value of the transverse response (nodal displacement in the z-direction) is evaluated after the insertion of 5 periodic cells, at a distance r = 6L from the cen- tre of the plate.As expected, few unit cells were enough to obtain a good vibration attenuation in the stop-band, even if the structure is finite.
The stop-band predicted using WFE method is located exactly at the place found using Floquet propagators with cylindrical wave correction; see Fig. 4b.In Fig. 6, this position is confirmed by a significant vibration attenuation in the stop-band detected by means of both the WFE and the standard FE models.The analysis of the response showed one of the main advantages of the WFE approach compared to standard FE approach: incredible reduction of the design and computational time due to the simplified structure and the significantly reduced size of the model.In this case study, for example, the size of the WFE model was a 12 × 12 matrix after the dynamic reduction, resulting in a WFE model of 6 DOFs, while the FE model was realised using 300364 triangular shell elements.
As seen in Fig. 5b, the location of the stop-band is weakly dependent upon the circumferential wavenumber m.Thus, the response to the time-harmonic excitation circumferentially distributed as cos(m ) , m = 1,2,3, remains qualitatively the same as for m = 0 (which is shown in Fig. 6) and is not presented.

Concluding Remarks
The reported results are summarised as follows: • The method of Floquet propagators with cylindrical wave correction and the WFE method are in a good agreement with each other in assessing location of stop-bands in radially periodic plates.Both methods, however, should be used for modelling of dynamics of a unit radial periodicity cell positioned at a distance, larger than its length, from the origin of polar coordinates.• Formulation of Floquet propagators with cylindrical wave correction for identification of stop-bands in radially periodic plates extends the theory of Bragg fibre to waveguides with constituents supporting both propagating and evanescent waves.Therefore, the method of Floquet propagators with cylindrical wave correction may be used for stop-band identification in a broad range of two-dimensional problems in polar and cylindrical coordinates, including high-order plate theories.
• Given the versatility of WFE method with respect to the library of available finite elements, its adjustment to deal with periodicity in polar coordinates facilitates its application for design and optimisation of modern acoustic metamaterials, featuring discrete and continuous constituents in spider net-shaped and similar arrays or lattices.
Advancing with coupled radial/circumferential periodicity in two-dimensional periodic systems and generalisation of both methods to study spherically periodic three-dimensional systems constitute the subjects of our ongoing work.

Appendix
whose solutions yield the relationship between the wavenumber and frequency (dispersion curves) and the displacement L (wave mode shapes) of the cross-section due to wave motion.

2D Formulation for Axisymmetric Structures
In this formulation, the projections k α and k r of the wavenumbers k in the radial and circumferential directions are considered.Imposing the periodicity conditions, the nodal DOFs = T 1 T on each side of the segment in Fig. 3 are related as where = e −ik and r = e −ik r L .Since the angle α is assumed to be small and the phase change of a wave as it propagates around the circumference must be a multiple of 2 , it must be Considering that is known for each circumferential order n, the problem is reduced to a 1D problem and the formulation presented in the first part of the appendix applies.

Fig. 1
Fig. 1 Schematic representation of the polar periodic plate model

Fig. 2
Fig. 2 Approximation of a periodic polar plate by a piecewise Cartesian waveguide and WFE discretisation of the unit cell j

Fig. 3
Fig. 3 Schematic representation of an infinite plate with periodically varying thickness, its unit cell, and WFE model