**1. Introduction**

Starting from the milestone book written by Leon Brillouin [1], wave propagation in periodic media has been a subject extensively studied. In particular, periodicity effects in electromagnetic wave motion have been largely investigated, and they have found applications in many optical and electromagnetic devices. The ability of periodic configuration of creating electronic/photonic band-gaps in semiconductors and crystals is similar to structural/acoustical band-gaps in elastic media. The subject has recently found renewed interest due to its new and potential applications in vibroacoustic isolation, e.g., [2], noise suppression devices, e.g., [3], mitigation of seismic waves, e.g., [4,5], elastic/acoustic metamaterials [6]. A comprehensive review of the research in periodic materials and structures has been presented in [7]. In [7] the authors gave an overview of the numerical and experimental research in periodic structures, phononic crystals, and acoustic/elastic metamaterials up to 2014, showing some of the recent progress and the growth in academic and applied research interests in these fields. Experimental studies on periodic structures were presented recently in many papers, and periodicity effects such as structural and acoustic band-gaps, attenuation, directional energy flow, were verified in many specific cases, e.g., [7,8]. While these phenomena are well known, most of the literature on periodic engi-

**Citation:** Manconi, E.; Sorokin, S.V.; Garziera, R.; Quartaroli, M.M. Free and Forced Wave Motion in a Two-Dimensional Plate with Radial Periodicity. *Appl. Sci.* **2021**, *11*, 10948. https://doi.org/10.3390/ app112210948

Academic Editor: Jin-Gyun Kim

Received: 21 September 2021 Accepted: 16 November 2021 Published: 19 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

neering structures has been devoted to developing theoretical and numerical approaches due to the importance of understanding and characterising wave propagation behaviour.

Mead [9] and his coworkers at the University of Southampton [10] gave a substantial contribution to methods for predicting and analysing wave motion in periodic engineering structures. Among the numerical approaches proposed, those based on Finite Element theory have shown large versatility and applicability to modelling physical structures. Two different main approaches are typically applied. The firsts implement absorbing boundary conditions/layers that mimick the infinity domains through artificial non-reflective boundaries [11,12]; the advantage of this approach relies on the capability of standard FEA to model waveguides as a part of more complex structures. However, the computational cost remains a common issue in such an analysis, especially at high frequencies, when the mesh must be refined to reach sufficient accuracy in studies of the performance of a waveguide. It implies a substantial increase in the number of DOFs of the model, which leads to prohibitive computational time at higher frequency.

The seconds are numerical methods based on FE analysis of a unit cell of the structures and the theory of wave propagation in periodic structures, e.g., [13–15]. Compared to FE modelling with absorbing boundary conditions/layers, these numerical methods allow high accuracy up to high frequency at very low computation cost, and they can be the preferable choice when the requirement is the prediction of wave motion in 1-dimensional and 2-dimensional single waveguides, e.g., beam-like structures, panels and shells, cylindrical waveguides, etc. Some FE commercial software have recently implemented new specific modules for studying Cartesian periodic structures exploiting Bloch-Floquet theory. Although their use can be advantageous due to their ability to tackle complicated geometries, they involve a high computational cost in many cases.

A common feature of these methods is Cartesian spatial periodicity. The Bloch-Floquet theorem, which was applied in 1946 by Brillouin to solve the wave equation, relies, in fact, on the translational invariance of the problem formulation and is not rigorously applicable to polar coordinates systems [16]. In many practical engineering situations, a source of vibrations may excite a large and flexible structure such as a ship's deck, an aeroplane fuselage, a satellite antenna, a wall panel. In these cases, radial periodicity (e.g., as a sequence of annuli with alternating properties) may be used to reduce transmission of the vibration and structure-borne sound, and numerical approaches to study free and forced wave propagation in polar coordinates are desirable. The problem of wave propagation in radially periodic structures has been mainly formulated in optics within the theory of Bragg fibre [17]. Leaving aside numerous purely numerical studies (see, for instance, [18–20]), we notice that two approaches have been used so far to "adjust" the Bloch theorem for a cylindrically symmetric Bragg fibre. The first one is based on the use of far-field approximation of Bessel functions [21,22], and the second one implies special radial varying of material properties [23–25]. Some recent works presented an approximation of the formulation of the Floquet theory for radially periodic membrane and plates [26]. However, studies of wave propagation in a polar periodic configuration in structural mechanics are rare.

The paper aims to present a numerical method based on Floquet-theory and the FE discretisation of a small finite slice of a radial periodic structure. The method relies on adapting the Wave Finite Element (WFE) approach [14,15] to radial periodic waveguides. In this method, a unit cell of the waveguide is discretised using standard FE elements. The FE mass and stiffness matrices are reduced using wave propagation theory in periodic structures, and wave characteristics are numerically evaluated from an eigenvalue problem. This allows a very substantial reduction in the number of DOFs involved in the computation, with a dramatic reduction of the computational time compared to other FE approaches. Forced wave propagation is then evaluated in the wave domain as described in [27,28]. Amongst the numerical methods that could be used to investigate wave propagation problems, the WFE technique has several desirable features: it can be applied both to continuous and periodic structures; it exploits Bloch-Floquet theory and the versatility of standard FE analysis of a very small part of the structure; it allows the study of waveguides

with complex cross-sectional characteristics in a systematic manner, up to high frequency and with a low-computation cost. Applications of the technique to periodic, axisymmetric, helical, and slowly varying waveguides were presented in the literature, and the method has been validated through many benchmark cases, e.g., [29–31]. The approach rigorously assumes wave propagation in structures posed in Cartesian coordinates, and application to forced response in two dimensions was only possible by a semi-analytical computation via contour integration and the residue theorem [32].

Following a previous study by the same authors [33], this paper presents a simplified adaptation of this WFE technique to structures in polar coordinates exhibiting periodicity in the radial directions. Cylindrical wave propagation is thus estimated exploiting the Floquet theory formulation for an infinite periodic structure in one dimension. The approximation is achieved by taking only a very small slice of the structure and discretising the slice through piecewise Cartesian segments. Wave characteristics in each segment are obtained by the WFE approach, while wave amplitudes change due to the changes in the geometry of the slice are accommodated in the model assuming that the energy flow through each segment is the same. In this paper, the forced response of a flexible periodic plate excited by a transverse source of vibration force is considered. Nearfield response and low-frequency harmonic excitation, far enough from the first stop-band, is assumed. The aim is to quantify the response level of the structure close to the excitation point and in the frequency regime where the fundamental modes of the structures are excited.

The paper is organised as follows. In Section 2 the method is described. Section 3 contains a numerical example of a thin flexible plate, for which analytical solutions are available [34]. Section 4 deals with a plate with a periodic change of thickness in the radial direction for which an analytical solution is not available. The numerical results were verified through a standard Finite Element model of the plate with perfectly matched layers (PMLs) in COMSOL Multiphysics®. The last section is devoted to conclusions.

#### **2. Stepwise WFE Approximation of a Radially Periodic Plate**

A lossless and linear elastic plate with radial periodicity is considered. The plate consists of an infinite repetition of a sequence of annuli of the same width with alternating properties, as shown in Figure 1a. The lattice constant, or the characteristic length of the unit cell, is here defined in the radial direction and denoted by *R*.

**Figure 1.** (**a**) radially periodic structure; (**b**) slice of the radially periodic structure; (**c**) approximation in piecewise Cartesian periodic waveguides.

In order to apply the WFE method, a slice of the plate is taken, as shown in Figure 1b. Since the periodicity is in the radial direction, the angle of the slice can be arbitrarily small. The slice is approximated as a piecewise rectangular waveguide in Cartesian coordinates, as shown in Figure 1c. A finite number of segments is assumed. These segments are numbered increasingly in the radial direction from 1 to *n*, where *n* is the cell distant *r* from the centre. Wave characteristics of each of these rectangular segments are obtained using the WFE approach as described in Section 2.1. In the model, the left nodes of the first cell are shifted at an arbitrarily small distance *R*<sup>0</sup> from the origin of the coordinates; the inner part, 0 < *r* < *R*0, is uniform and homogenous.

#### *2.1. Free Wave Propagation Characteristics*

Wave characteristics of each segment (cell) in Figure 1c are evaluated using the WFE method as described in [14,28]. The waves are numerically represented by the dispersion curves, (**k**, *ω*), which give the information about the wave vector **k** available for each frequency *ω*, and the corresponding FE nodal displacements **Φ<sup>q</sup>** and nodal forces **Φf**, which occur under the passage of a wave. In the WFE, only the unit periodic cell is analysed. Figure 2 shows a finite element mesh of a unit cell using 4-noded rectangular elements and the correspondent WFE model. The degrees of freedom **q** in the dicrestised equation of motion, Equation (1), are ordered and condensed as **q** = # **q***T <sup>L</sup>* **<sup>q</sup>***<sup>T</sup> <sup>R</sup>* **<sup>q</sup>***<sup>T</sup> I* \$*T* , where the subscripts *L, R* and *I* are associated with the right **q***<sup>R</sup>* = # **q***T* <sup>1</sup> **<sup>q</sup>***<sup>T</sup>* 2 \$*T* , left **q***<sup>L</sup>* = # **q***T* <sup>5</sup> **<sup>q</sup>***<sup>T</sup>* 6 \$*T* and interior **q***<sup>I</sup>* = # **q***T* <sup>3</sup> **<sup>q</sup>***<sup>T</sup>* 4 \$*<sup>T</sup>* nodal degrees of freedom of the cell, with a similar expression for the nodal forces **f** = ! **f** *T <sup>L</sup>* **f** *T <sup>R</sup>* **f** *T I* "*T* .

**Figure 2.** (**a**) Finite element mesh of a unit cell using 4-noded rectangular elements and (**b**) WFE nodes condensation.

The discretised equation of motion, (**<sup>K</sup>** − *<sup>ω</sup>*2**M**)**<sup>q</sup>** = **Dq** = **<sup>f</sup>**, is obtained using conventional FEA and it is partitioned as

$$
\begin{bmatrix}
\mathbf{D}\_{LL} & \mathbf{D}\_{LR} & \mathbf{D}\_{RI} \\
\mathbf{D}\_{RL} & \mathbf{D}\_{RR} & \mathbf{D}\_{RI} \\
\mathbf{D}\_{IL} & \mathbf{D}\_{IR} & \mathbf{D}\_{II}
\end{bmatrix}
\begin{Bmatrix}
\mathbf{q}\_{L} \\
\mathbf{q}\_{R} \\
\mathbf{q}\_{I}
\end{Bmatrix} = \begin{Bmatrix}
\mathbf{f}\_{L} \\
\mathbf{f}\_{R} \\
\mathbf{f}\_{I}
\end{Bmatrix} \tag{1}
$$

To apply the WFE method, the internal Dofs of the unit cell must be reduced. In this paper, a standard dynamic condensation is applied as

$$\mathbf{q}\_I = -\mathbf{D}\_{II}{}^{-1}(\mathbf{D}\_{IL}\mathbf{q}\_L + \mathbf{D}\_{IR}\mathbf{q}\_R) \tag{2}$$

However, according to the size of the model, more efficient methods for condensing the inner DOFs and speeding up the resolution of the WFE eigenvalue problem can be applied, e.g., [35,36].

The resulting equation of motion is

$$
\begin{bmatrix}
\mathbf{D}\_{LL} & \mathbf{D}\_{LR} \\
\mathbf{D}\_{RL} & \mathbf{D}\_{RR}
\end{bmatrix}
\begin{Bmatrix}
\mathbf{q}\_L \\
\mathbf{q}\_R
\end{Bmatrix} = \begin{Bmatrix}
\mathbf{f}\_L \\
\mathbf{f}\_R
\end{Bmatrix} \tag{3}
$$

The inclusion of stress, temperature effects and damping in the model can be accommodated, if necessary, as described in [37]. Periodicity conditions are applied so that

$$\mathbf{q}\_R = \lambda \mathbf{q}\_L; \quad \mathbf{f}\_R = -\lambda \mathbf{f}\_L \tag{4}$$

where *λ* = *e*−*ikL* and *k* is the wavenumber. Combining Equation (4) with Equation (3), the equation of free wave motion takes the form of a quadratic eigenvalue problem, as

$$
\left[\lambda^2 \mathbf{D}\_{\rm LR} + \lambda (\mathbf{D}\_{\rm LL} + \mathbf{D}\_{\rm RR}) + \mathbf{D}\_{\rm RL}\right] \mathbf{q}\_{\rm L} = 0 \tag{5}
$$

whose solutions yield the relationship between the wavenumber and frequency (dispersion curves) and the displacement **q***<sup>L</sup>* of the cross-section due to wave motion (wave mode shapes). Equation (5) can be further recast as a standard linear eigenvalue problem as (**A** − *γ***I**)**Z** = 0,

$$\begin{aligned} \text{where } \mathbf{A} &= \begin{bmatrix} 0 & \mathbf{I} \\ -\mathbf{A}\_2^{-1} \mathbf{A}\_0 & -\mathbf{A}\_2^{-1} \mathbf{A}\_1 \end{bmatrix} \mathbf{A}\_2 = \mathbf{D}\_{LR}, \mathbf{A}\_1 = \mathbf{D}\_{LL} + \mathbf{D}\_{RR}, \mathbf{A}\_0 = \mathbf{D}\_{RL}, \text{ and} \\ \mathbf{Z} &= \begin{bmatrix} \underline{\mathbf{q}}\_L^T & \lambda \mathbf{q}\_L^T \end{bmatrix}^T. \end{aligned}$$

The time average energy flow in each segment can be then obtained from

$$
\Pi = \frac{1}{2} \text{Re} \left[ \mathbf{f}^H i\omega \mathbf{q} \right] \tag{6}
$$

where the superscript *H* denotes the Hermitian transpose and **f** and **q** are recovered using Equations (3) and (4).

#### *2.2. Forced Wave Amplitude*

Forced wave propagation and structural response are evaluated following the theory presented in [27,28]. Wave properties are grouped into positive and negative going waves, which can be described by the two sets of results **k**+, **Φ**+, **a**+ and **k**−, **Φ**−, **a**− . Here **k+**, **k**<sup>−</sup> and **a+**, **a**<sup>−</sup> are wavenumbers and waves' amplitudes travelling in the positive and negative direction, while **Φ**<sup>+</sup> = [**Φ**<sup>+</sup> **<sup>q</sup>** , **<sup>Φ</sup>**<sup>+</sup> **f** ] *<sup>T</sup>* and **Φ**<sup>−</sup> = [**Φ**<sup>−</sup> **<sup>q</sup>** , **Φ**<sup>−</sup> **f** ] *<sup>T</sup>* are the corresponding nodal displacements and nodal forces, that is the FE discretisation of wavemodes. For evaluating the forced response, it is advantageous to obtain also the left eigenvectors of the WFE eigenvalue problem in Equation (5). These are partitioned in the same manner as the wavemodes, that is **Ψ**<sup>+</sup> = [**Ψ**<sup>+</sup> **<sup>f</sup>** , **<sup>Ψ</sup>**<sup>+</sup> **<sup>q</sup>** ] and **Ψ**<sup>−</sup> = [**Ψ**<sup>−</sup> **<sup>f</sup>** , **Ψ**<sup>−</sup> **<sup>q</sup>** ]. Left and right eigenvectors are orthogonal and they can be normalised so that

$$\mathbf{Y}^+\Phi^+ = \mathbf{I}, \mathbf{Y}^-\Phi^- = \mathbf{I} \tag{7}$$

In an analogous manner to the modal analysis, the total displacement and force at the junction of a cell is described by the sum of the positive and negative wavemodes so that

$$
\begin{bmatrix} \mathbf{q} \\ \mathbf{f} \end{bmatrix} = \begin{bmatrix} \Phi\_{\mathbf{q}}^{+} & \Phi\_{\mathbf{q}}^{-} \\ \Phi\_{\mathbf{f}}^{+} & \Phi\_{\mathbf{f}}^{-} \end{bmatrix} \begin{bmatrix} \mathbf{a}^{+} \\ \mathbf{a}^{-} \end{bmatrix}. \tag{8}
$$

Equation (8) defines the transformation between the physical domain, where the motion is described in terms of nodal displacements and forces, and the wave domain, where the motion is described in terms of waves of amplitudes **a+**, **a**<sup>−</sup> travelling in the positive and negative directions. A point force **f***e* will generate excited positive and negative going waves, propagating away from the excitation point. Compared to the Cartesian periodicity, in the case of a radially periodic plate, the position of this point force is an independent parameter. Here, we assume a transverse harmonic point force exciting the structure at the left nodes of the first cell in Figure 1c. With reference to Figure 1, continuity

and equilibrium equations can be written at the left side of the first cell, and excited wave amplitudes can be recovered from Equation (9)

$$
\begin{bmatrix}
\Phi\_{\mathbf{q},1}^{+} & \Phi\_{\mathbf{q},1}^{-} \\
\Phi\_{\mathbf{f},1}^{+} & \Phi\_{\mathbf{f},1}^{-}
\end{bmatrix}
\begin{bmatrix}
\mathbf{a}\_{1}^{+}(r\_{1}) \\
\mathbf{a}\_{1}^{-}(r\_{1})
\end{bmatrix} = \begin{bmatrix}
0 \\
\mathbf{f}\_{\mathbf{f}}
\end{bmatrix} \tag{9}
$$

In practice, as in modal analysis, only *m* pairs of (positive- and negative-going) waves are retained (a reduced wave basis is assumed in most cases). Moreover, the number of assumed modes can be different at each frequency. It is noteworthy to mention that all the waves propagating in the structure (real wavenumbers) should be considered. Evanescent and attenuating waves (corresponding to pure imaginary or complex wavenumbers) should be also assumed since they play an important role in scattering and forcing problems. Consequently, the wavemodes' matrix is typically rectangular and the use of standard pseudoinverse operation can lead to numerical errors. To overcome this issue, the orthogonality between right and left eigenvectors in Equation (7) can be exploited, and the excited wave amplitudes at *r*<sup>1</sup> can be obtained from

$$
\begin{bmatrix}
\mathbf{a}\_1^+(r\_1) \\
\mathbf{a}\_1^-(r\_1)
\end{bmatrix} = \begin{bmatrix}
\mathbf{Y}\_{f,1}^+ & \mathbf{Y}\_{\mathbf{q},1}^+ \\
\mathbf{Y}\_{\mathbf{f},1}^- & \mathbf{Y}\_{q,1}^-
\end{bmatrix} \begin{bmatrix}
0 \\
\mathbf{f}\_{\mathbf{f}}
\end{bmatrix} \tag{10}
$$

## *2.3. Coupling of the Segments and Wave Amplitude Decay*

Compared to the corresponding Cartesian waveguides, wave propagation is not translational invariant in the radially periodic structures, and amplitude attenuation occurs as the waves travel from the centre in the radial direction. Since the energy flowing along the slice is constant, wave amplitude change due to changes in the geometry can be accommodated in the model according to an energy balance principle as in [31]. In this section, the procedure is presented for the first two cells. The same passages must be applied up to the cell *j* = *n*, including the point at which the response must be evaluated.

In the Cartesian cell, amplitudes are related at two points *r*<sup>1</sup> and *r*<sup>2</sup> = *r*<sup>1</sup> + *R*, by **a**1 +(*r*2) = **T**+(*R*)**a**<sup>+</sup> <sup>1</sup> (*r*<sup>1</sup> and **a**− <sup>1</sup> (*r*2) = **T**<sup>−</sup> <sup>1</sup> (*R*)**a**<sup>−</sup> <sup>1</sup> (*r*<sup>1</sup> , where **T**<sup>+</sup> <sup>1</sup> (*R*) = diag# exp <sup>−</sup>i**k**<sup>+</sup> <sup>1</sup> *R* \$ and **T**− <sup>1</sup> (*R*) = diag# exp −i**k**<sup>−</sup> <sup>1</sup> *R* \$. Therefore, nodal displacements and nodal forces at the interface between cell 1 and cell 2 are

$$
\begin{bmatrix} \mathbf{q}\_1 \\ \mathbf{f}\_1 \end{bmatrix} = \begin{bmatrix} \Phi\_{\mathbf{q},1}^+ & \Phi\_{\mathbf{q},1}^- \\ \Phi\_{\mathbf{f},1}^+ & \Phi\_{\mathbf{f},1}^- \end{bmatrix} \begin{bmatrix} \mathbf{T}\_1^+ (\overline{\mathcal{R}}) & 0 \\ 0 & \mathbf{T}\_1^- (\overline{\mathcal{R}}) \end{bmatrix} \begin{bmatrix} \mathbf{a}\_1^+ (r\_1) \\ \mathbf{a}\_1^- (r\_1) \end{bmatrix} \tag{11}
$$

Using the left eigenvectors as in Equation (10), the first attempt to find the excited wave amplitude in cell 2 gives

$$
\begin{bmatrix} \stackrel{\sim}{\mathbf{a}\_2^+} \\ \stackrel{\sim}{\mathbf{a}\_2^-} \end{bmatrix} = \begin{bmatrix} \mathbf{Y}\_{f,2}^+ & \mathbf{Y}\_{q,2}^+ \\ \mathbf{Y}\_{f,2}^- & \mathbf{Y}\_{q,2}^- \end{bmatrix} \begin{bmatrix} \mathbf{q}\_1 \\ \mathbf{f}\_1 \end{bmatrix} \tag{12}
$$

Nodal displacements and nodal forces at cell 2 can be obtained using the same expression of Equation (10)

$$
\begin{bmatrix} \mathbf{q}\_2 \\ \mathbf{f}\_2 \end{bmatrix} = \begin{bmatrix} \Phi\_{\mathbf{q},2}^+ & \Phi\_{\mathbf{q},2}^- \\ \Phi\_{\mathbf{f},2}^+ & \Phi\_{\mathbf{f},2}^- \end{bmatrix} \begin{bmatrix} \mathbf{T}\_2^+(\overline{\mathcal{R}}) & 0 \\ 0 & \mathbf{T}\_2^-(\overline{\mathcal{R}}) \end{bmatrix} \begin{bmatrix} \tilde{\mathbf{a}}\_2^+ \\ \tilde{\mathbf{a}}\_2^- \end{bmatrix} \tag{13}
$$

and the time-averaged energy flows, Equation (6), of cell 1 and cell 2 can be evaluated using nodal displacement and nodal forces from Equations (11) and (13)

$$\Pi\_1 = \frac{1}{2} \text{Re} \left\{ \mathbf{f}\_1^H i\omega \mathbf{q}\_1 \right\}, \quad \Pi\_2 = \frac{1}{2} \text{Re} \left\{ \mathbf{f}\_2^H i\omega \mathbf{q}\_2 \right\} \tag{14}$$

The ratio between Π<sup>1</sup> and Π<sup>2</sup> gives the amplitude decay *ξ*<sup>2</sup> in waveguide 2, which is evaluated as in [32]

$$\zeta\_2^{\tau\_2^{\tau\_2}} = \frac{\text{Re}\{\mathbf{f}\_1^H \dot{i}\omega \mathbf{q}\_1\}}{\text{Re}\{\mathbf{f}\_2^H \dot{i}\omega \mathbf{q}\_2\}}\tag{15}$$

The wave and amplitudes in waveguide 2 are therefore approximated to

$$
\begin{bmatrix}
\mathbf{a\_2}^\* \\
\mathbf{a\_2}^- \\
\end{bmatrix} = \tilde{\zeta}\_2 \begin{bmatrix}
\tilde{\mathbf{a\_2}}^\* \\
\tilde{\mathbf{a\_2}}^- \\
\end{bmatrix} \tag{16}
$$

Equations (13)–(16) can be repeated iteratively until *ξ*<sup>2</sup> converges to one, and the final value of the wave amplitudes reach the required approximation. One passage was sufficient to converge to a useful approximation with a very low computational cost in the numerical cases studied. Excited wave amplitudes decay in the next cells up to cell *n* are evaluated following the same passages.

#### **3. Numerical Examples**

This section shows two numerical examples: the first concerns a thin isotropic plate excited by a central transverse harmonic force, while the second deals with a plate with a periodic radial change of thickness. In both cases, the WFE model is obtained using 4-noded plane elements in bending having three degrees of freedom per node: translation in the *z*-direction and rotations around the *x*- and *y*-directions. To verify the results, analytical solutions were compared in the first case, while results for the radially periodic plate were verified through comparison with those obtained by an FE model with PMLs. This model was realised in COMSOL Mutiphysics® using the Structural Mechanics Module. The plate was discretised by shell elements having six degrees of freedom per node (translations and rotations in the *z, x*, and *y* directions), while a second-order polynomial stretching function was chosen for the PMLs, see the COMSOL Multiphysics Reference Manual for further information. The FE models were realised both in Cartesian and Polar coordinates systems with similar results, and a convergence test was performed by refining the mesh in the PML domain.

The non-dimensional frequency

$$
\Omega = \omega \overline{R}^2 \sqrt{\frac{\rho h}{D}} \tag{17}
$$

is introduced, where *ρ*, *h* and *D* are respectively the density, thickness and bending stiffness of the plate, *ω* is the frequency in radiant, and *R* is the period of the structure as in Figure 1.

#### *3.1. Numerical Verification. Infinite Plate Subjected to a Transverse Harmonic Force*

The literature has largely studied the problem of a thin plate subjected to harmonic loads, and closed-form solutions can be found in many classical books on vibrations, e.g., [34]. This numerical example is introduced here to verify the method described in Section 2 and show its applicability to continuous structures.

An aluminium plate of thickness *h* = 1 mm is considered. Flexural waves in the plate are excited by a central transverse harmonic force of magnitude *F*/ *Eh*<sup>2</sup> = 1.4 · <sup>10</sup><sup>−</sup>3, and the complex out of plane displacement is evaluated at a distance *r*/*h* = 300 from the plate centre, viz. origin of the polar coordinates.

In this case, the structure is continuous and uniform, and it can be studied as a periodic structure with arbitrary radial and circumferential periodicity: the choice of the period, viz. length *R* of the unit cell, is arbitrary under standard FE assumptions to avoid element distortion and dispersion errors [38]. Therefore, the number of segments that are used to approximate the slice from the centre to *r* is arbitrarily chosen. The discretisation must be refined until it does not produce a negligible change in the solution. In this example, only three segments have been found sufficient for convergence.

In order to simplify the WFE model, a small circle around the point force is removed and the left nodes of the first cell are shifted at a distance *R*0/*h* = *r*<sup>1</sup> = 10 from the point force. In this example, the WFE approximation is obtained for a slice of angle *θ* = 6o. Waves are induced only in the positive radial direction and therefore **a**<sup>1</sup> − = 0 is assumed. Periodicity in each segment is further exploited. Here, the number of "sub-period" in each segment is chosen by a simple algorithm that optimises the mesh according to the slice's dimension and the number of segments. No significant differences are noticed decreasing the distance *R*<sup>0</sup> further.

Figure 3 shows the absolute value and the phase of the transverse displacement of the plate. A comparison between the numerical and analytical solutions is shown. The latter is evaluated using the analytical equation given in [34]: *w* = <sup>−</sup>*iF* 8*Dk*<sup>2</sup> *H*<sup>2</sup> <sup>0</sup> (*kr*) <sup>−</sup> *<sup>i</sup>* <sup>2</sup> *<sup>π</sup> K*0(*kr*) , where *k* = \*<sup>4</sup> *ω*2*ρh*/*D* is the flexural wavenumber, *H*<sup>2</sup> <sup>0</sup> (*kr*) is the zero-order Hankel function of the second kind and *K*0(*kr*) is the zero-order modified Bessel function of the second kind.

**Figure 3.** Response of a thin plate excited by a harmonic transverse point force, (**a**) absolute value; (**b**) phase. Comparison between the approximated WFE, FE model with PMLs and analytical results.

Results obtained using the FE model of the plate with PMLs are, as seen in Figure 3, in excellent agreement with the WFE results (small discrepancies are due to the differences in the FE models). However, this numerical study also demonstrates the computational superiority of the WFE method over the standard FE + PML analysis: the CPU time used by Matlab® to solve the WFE model of the plate was approximately four hundred times less than the CPU time used by COMSOL Mutiphysics® to solve the correspondent FE + PML model. One of the main reasons for such a difference in the computational cost is the very small size of the WFE matrices compared to the number of DOFs involved in standard FE analysis. The WFE model for this numerical example was realised using one shell element, resulting in 6 DOFs after the WFE reduction, Equation (5), while the size of the FE + PML model was 676,422 DOFs.

#### *3.2. Radially Periodic Plate*

In this section, we consider a plate in polar coordinates with radial periodicity. Figure 4 depicts a schematic figure of the plate and its unit cell (period). In the following, the characteristic length *R* and the thickness *h*<sup>1</sup> are used to define the dimensionless parameters:

$$\beta = h\_2/h\_1 = 2.7, \ \gamma = \overline{\text{R}}/R\_1 = 2.7, \ \delta = \overline{\text{R}}/h\_1 = 27$$

**Figure 4.** (**a**) radially periodic infinite plate with an internal hole of radius *R0*; (**b**) period, or unit cell.

The WFE model is obtained considering an arbitrary circumferential periodicity of angle *θ* = 9o. The external harmonic force of magnitude *F*/(2*πR*0) =140 N/m is applied at the internal edge of the plate and included in the WFE by two equivalent transverse nodal forces exciting the left nodes of the first unit cell. Waves are induced only in the positive radial direction and therefore **a**<sup>1</sup> − = 0 is assumed. The structure exhibits stop-bands due to its periodicity as shown in [26]. These stop-bands can be predicted with accuracy by the WFE up to high frequency. However, in this paper, the plate response is evaluated in the nearfield and at a low-frequency. Therefore, only the wave modes below the first stop-band are considered in the analysis. Figure 5 shows the dispersion curves of flexural waves propagating in the positive direction. Analytical flexural waves propagating in the corresponding homogeneous plates of thickness *h*<sup>1</sup> and *h*<sup>2</sup> are also shown for comparison.

**Figure 5.** Dispersion curves for flexural waves in the radially periodic plate. Propagating flexural waves in the corresponding homogeneous plates of thickness *h*1, **----** , and *h*2, **-.-.-** , are shown for comparison.

The transverse response of the plate is evaluated after 4 unit cells, at a distance *r* = 5*R* from the centre of the plate (or equally at a distance *r* = 4*R* from the border of the internal hole). Figure 6 shows the displacements in terms of absolute and phase values. Results obtained using the FE model with PMLs are also presented. Although some small discrepancies were expected due to the very different FE discretisation and a very different approach to the problem, it can be noticed that the results obtained by the WFE method are in good agreement with those obtained by the finite element harmonic analysis. The WFE model for each segment of the plate slice was set up with three shell elements, resulting in six DOFs after the WFE reduction, while the FE + PML had 794,448 DOFs. The

computational time used by Matlab® to solve the WFE model of this radially periodic plate was almost two thousand times less than the CPU time used by COMSOL Mutiphysics® to solve the correspondent FE + PML model.

**Figure 6.** Response at the 4th cell of the radial periodic plate excited by harmonic transverse point forces at the border of the central hole, (**a**) absolute value; (**b**) phase. Comparison between the approximated WFE results FEM results.

#### **4. Conclusions**

In this paper, free and forced wave propagation in a radially periodic plate was studied using an adaptation of the Wave Finite Element method to structures showing radial periodicities. Cylindrical waves propagating were approximated exploiting the Floquet theory formulation for an infinite periodic structure in one dimension. This approximation was achieved by taking only a very small slice of the structure, which was discretised through piecewise Cartesian segments. Wave characteristics in each segment were obtained by the WFE method, while wave amplitudes change was accommodated in the model assuming an energy balance principle. The forced response of the structure was then evaluated in the wave domain. The paper's main goal was to verify the approach for predicting the forced response in the nearfield of a radial periodic structures. Two numerical examples were presented: an isotropic thin plate excited by a transverse harmonic force, for which analytical solutions are available, and a plate with a periodic change of thickness in the radial direction. In the latter, the numerical results were verified through a standard Finite Element model of the plate with perfectly matched layers (PMLs). In both cases, it was found that the numerical results were in good agreement with the analytical and numerical FE results, showing the advantages of the approach in terms of computational ime, approximation controllability, and modelling efficiency.

**Author Contributions:** Conceptualization, E.M., S.V.S. and R.G.; methodology, E.M. and S.V.S.; software, E.M. and M.M.Q.; validation, R.G. and M.M.Q.; formal analysis, E.M.; investigation, E.M. and M.M.Q.; data curation, E.M.; writing—original draft preparation, E.M. and S.V.S.; writing review and editing, E.M. and S.V.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

