## Abstract

This paper presents a mathematical model describing propagation of bending waves in a perforated thin plate. It is assumed that the holes are circular and form a doubly periodic square array. A spectral problem for the biharmonic operator is formulated in a unit cell containing a single defect, and its analytical solution is constructed using a multipole method. The overall system for the coefficients in the multipole expansion is then solved numerically. We generate dispersion diagrams for the two cases where the boundaries of holes are either clamped or free. We show that in the clamped case, there is a total low-frequency band gap in the limit of inclusions of zero radius, and give a simple formula describing the corresponding band diagram in this limit. We show that in the free-edge case, the band diagram of the vibrating plate is much closer to that of plane waves in a uniform plate than for the clamped case.

## 1. Introduction

This paper addresses a new analytical model of flexural Bloch–Floquet waves in an infinite plate containing a doubly periodic array of circular holes (with clamped or free edges).

Waves in homogeneous plates are well described in the classical literature (e.g. Graff 1975), while the analysis of flexural waves in heterogeneous periodic structures was initiated very recently. A recent paper by Halkjær *et al.* (2006) addresses the problems of numerical optimal design and maximizing the band gaps for periodic plate structures. The analytical work of Evans & Porter (2007) includes the analysis of flexural waves in thin plates periodically constrained at a square array of points of an infinite plane.

In the present paper, we construct a full analytical solution of the spectral problem for the biharmonic operator in an infinite plane with a square array of circular holes. Both ‘free-edge’ and ‘clamped-edge’ boundary conditions are addressed here. The multipole expansions, written in terms of the Bessel and modified Bessel functions, are used to represent the eigensolutions. The dispersion equations are constructed in an explicit analytical form and the results of the accompanying numerical simulations are compared with independent finite-element simulations.

Among other results, we demonstrate that the clamped system of circular cylinders possesses a full band gap at zero frequency, and that this band gap persists even when the radius of the cylinders tends to zero. This means that the Dirichlet problem for flexural waves in two dimensions cannot be homogenized. On the other hand, we show that holes with free boundaries change the properties of flexural waves only in a perturbative fashion.

The structure of the paper can be described as follows. Section 2 gives the mathematical formulation of the flexural wave problem on an elementary cell. The matrix representation of the boundary conditions, in terms of the multipole coefficients, is discussed in §3. Rayleigh identities and lattice sums are derived in §4. The derivation of the algebraic formulation for the eigenvalue problem is completed in §5. Section 6 presents the results of analysis of the dispersion equation and numerical simulations of eigenfields for the clamped- and free-edge boundary conditions.

## 2. Formulation of the spectral problem for a biharmonic operator in a unit cell

We consider a doubly periodic square array of circular voids within a thin elastic plate, as shown in figure 1. The elementary central cell of the periodic structure is defined bywhere

The out-of-plane elastic displacement satisfies the equation of motion(2.1)where *t* is the time; *h* is the plate thickness; *ρ* is the mass density; is the flexural rigidity of the plate; *E* denotes the Young modulus; and *ν* is the Poisson ratio of the elastic material.

We are interested in the time-harmonic vibrations of the plate, so thatwhere (*r*, *θ*) are the polar coordinates and *p* is the angular frequency. In this case, equation (2.1) is transformed to the form(2.2)where

The solution of (2.2) admits the representation(2.3)where(2.4)and hence *W* is the superposition of ‘waves’ with real and imaginary wave numbers. In the text below, we shall use the series representation(2.5)with *W*_{n} defined by(2.6)where, in (2.6), *J*_{n} is the Bessel function; is the Hankel function; and *I*_{n} and *K*_{n} are the modified Bessel functions related to *J*_{n} and , respectively, by (see Abramowitz & Stegun 1965, eqns 9.6.3 and 9.6.4)(2.7)Note that the expansion (2.5) converges only in the region between (say) the cylinder at the origin and the circle touching its four nearest neighbours: . This means that the expansion covers the entire unit cell if . For radii larger than this, note first that our main result from which band diagrams are calculated is the Rayleigh identity, which requires for its derivation only that there be a circle lying between the boundaries of both the void and the unit cell (i.e. *r*<*d*/2). Second, if it is necessary to calculate properties associated with the solution throughout the unit cell for such larger radii, then the single origin expansion (2.5) needs to be replaced by a multiple origin expansion incorporating only Bessel functions of the types and . Such an expansion has been referred to as a Wijngaard expansion (Chang & Chang 1994; Lo *et al.* 1994), but it is not required in the numerical examples here.

The Bloch–Floquet quasi-periodicity conditions are set within the periodic array, and they relate the solution within the general cell to that in the central cell as(2.8)where(2.9)*k*_{0} is the Bloch vector; *e*^{(1)} and *e*^{(2)} are the axial unit vectors; and *d* is the period of the structure.

We will consider two types of boundary conditions on the circular contour , which are as follows.

Clamped boundary(2.10)

Free boundary where the moment and the transverse force are both equal to zero, i.e.(2.11)where(2.12)

(2.13)and(2.14)

Our aim is to solve the above spectral problem by constructing a dispersion relation, which links the angular frequency *p* with the Bloch vector *k*_{0}. The coefficients in the multipole representations (2.5) and (2.6) will be evaluated as solutions of the Rayleigh system, which will be derived in §4.

## 3. Multipole representation of the boundary conditions

We first consider the simpler case of the clamped boundary. Equations (2.6) and (2.10) imply(3.1)and(3.2)

Next, we address the case of the free boundary. The first condition in (2.11) together with (2.6) gives(3.3)which (with the use of the standard Bessel equations for and , respectively) can be transformed to the form(3.4)The second condition in (2.11) takes the form(3.5)or equivalently,(3.6)Equation (3.6) can be simplified using standard Bessel equations, which give(3.7)and(3.8)and similarly for the other two types of Bessel functions. Thus, we obtain(3.9)

## 4. Rayleigh identities and lattice sums

According to (2.3) and (2.4), the functions *W*^{(1)} and *W*^{(2)} are oscillatory and exponential, respectively, and propagate independently in the elastic medium between the periodic circular voids (inclusions). Hence, they have separate consistency relations between the local expansions in the elastic medium surrounding a general inclusion and the inclusion in the central cell. These consistency relations are derived in the standard way (e.g. Botten *et al.* 2005) using Green's theorem and Graf's addition theorem (Abramowitz & Stegun 1965), and are called the Rayleigh identities. Note that Abramowitz and Stegun give only the addition theorem for Bessel functions of the types *J*, *Y*, *H*^{(1)} and *H*^{(2)}, collectively denoted by the symbol . Their formula (9.1.79) is(4.1)where the vectors *u*, *v* and *w* lie on the edges of a triangle, with *α* the angle between *u* and *v*, and *Χ* the angle between *u* and *w*, *v* being shorter than *u* . We will also need the addition theorem for the modified Bessel functions *K* and *I*. This is obtained by replacing *u*,*v*,*w* by , and using the connection between and , and between and . The result is(4.2)with a similar formula for *I*_{n}.

The Rayleigh identity for the function *W*^{(1)}, which satisfies the Helmholtz equation (2.4), is given by(4.3)In (4.3), the quantity is the lattice sum of order *l* for the Hankel function *H*^{(1)}(4.4)where *R*_{q} is defined in (2.9) and has the polar coordinate representation Note that the sums in (4.4) are not absolutely convergent, so that, while the representation is useful analytically, it needs to be replaced by alternative expressions for numerical applications.

We also have (Chin *et al.* 1994), where various ways of evaluating the non-trivial lattice sums are reviewed by Moroz (2006). Here, we calculate using the representation obtained by repeated integration to accelerate convergence (Chin *et al.* 1994)(4.5)Here, denotes a pair of integers, which are summed over a sufficiently large range to guarantee the desired accuracy, while runs over vectors with the required quasi-periodicity of the Bloch vector *k*_{0}. The parameter *ξ* is arbitrary and is generally set equal to *d*, while the acceleration parameter *p* is chosen in the range 3–7 typically (increasing it makes the series converge more rapidly for large *Q*_{h}, but delays the onset of this rapid convergence).

The Rayleigh identity for the function *W*^{(2)}, which satisfies the modified Helmholtz equation (second equation in (2.4)), has the form(4.6)where(4.7)which is exponentially convergent and hence can be evaluated by direct summation for *βd* not very small. The non-accelerated version of these lattice sums corresponding to (4.5) is(4.8)The version of this corresponding to acceleration by threefold integration is(4.9)This is applied for the region of small *β*, which will be important later for the case of holes with free edges.

Note that accelerated representations like (4.5) and (4.9) can only be applied for *l*>0. To use them for *l*<0, we need relations between the sums of positive and negative orders. The following identities hold for any lattice:(4.10)and(4.11)where the asterisk denotes complex conjugation.

## 5. Homogeneous matrix equation for the eigenmodes of elastic vibrations

The overall system of equations combines the Rayleigh identities (4.3) and (4.6) and the boundary conditions (3.1), (3.2) or (3.4), (3.8).

We note that the Rayleigh identities couple multipoles of different orders for the same type of waves and depend on the geometrical arrangement of the inclusions within the periodic array. By contrast, the boundary conditions (3.1), (3.2) and (3.4), (3.8) do not couple different multipole orders, but couple multipoles of the same order for the two types of waves. They contain the information about the elastic moduli, the radius of the inclusion and the type of the boundary condition.

Let *L* denote the order of truncation of the multipole series (2.5), so that the multipole summation runs from −*L* to *L*. Thus, the lattice sums must be calculated over orders ranging from to

Consider first the case of a clamped boundary. Substitution of (4.3) and (4.5) into (3.1) and (3.2) gives(5.1)and(5.2)

For the truncated system, arranging the unknown coefficients and into column vectors ** E** and

**, equations (5.1) and (5.2) can be combined into a single matrix equation(5.3)where each of the four sub-blocks has the dimension**

*F*Next, we consider the case of a free boundary. Here, we eliminate by substitution of (4.3) and (4.6) into the boundary conditions (3.4) and (3.6). This gives(5.4)and(5.5)As in the previous case, the matrix form of the truncated system (5.4) and (5.5) is(5.6)where each of the sub-blocks again has the dimension

The equations for holes with clamped (5.3) or free (5.6) edges take the formThe dispersion equation then iswhich we solve numerically for a given value of *k*_{0} for the desired number of solutions *β*. The truncation order *L* we use here is 1, so the complex-valued matrix has dimension 6×6. Larger values of *L* can be used when the holes are closer to touching each other than in the cases studied here.

## 6. Asymptotic approximations and band diagrams

We commence with the case of a square array (with a period unity) of holes whose boundaries are clamped. We give a band diagram for this problem in figure 2. This diagram follows the usual convention from the photonic band gap literature, in that three principal lines are chosen to represent the variation of : from to ; from to ; and (to the left) from to , with . The distance from the origin along the horizontal axis is the modulus of the Bloch vector *k*_{0}. Features to be noticed in figure 2 are the total band gaps above and below the fundamental mode, and the absence of a band gap between the second and third modes.

A very interesting feature of the diagram in figure 2 is the band gap at low frequencies. A similar feature occurs in electromagnetism for arrays of perfectly conducting rods (Nicorovici *et al.* 1995*a*,*b*; Pendry *et al.* 1996; Guida *et al.* 1998), where it is important in providing one of the essential elements required for the construction of metamaterials.

We can obtain an expression giving an estimate for the *β*-value of the first band near normal incidence by truncating the system (5.1) and (5.2) to include only terms .

The determinant of the resulting 2×2 matrix is(6.1)For the data of figure 2, with and , the zero of this determinant occurs at a *β*-value of 5.638, compared with the zero of the determinant of the 6×6 matrix, which occurs at 5.644. The multipole method used here converges very rapidly for cylinders of small radius, as evidenced by this example.

In figure 3, we compare the result of the multipole method (with ) and an independent finite-element computation (Comsol v. 3.2) for the displacement *W* of a thin plate perforated by circular holes with clamped boundaries, with the Bloch factor prescribed at and . Note that the mode scale factor is arbitrary and is different for the two calculations, but both give very similar displacement maps. The mode here is symmetric and well described by a monopole approximation , as just noted. (Small numerical errors in *W* near the corners of the unit cell are apparent in figure 2, in keeping with the comments on convergence in §2.)

We can take the limit as the radius *a* tends to zero in (6.1). The result is the following simple equation, which gives the band diagram of a square lattice of clamped inclusions of zero radius:(6.2)The band diagram corresponding to (6.2) is given in figure 4, together with lines marking the singularities of the lattice sum . Note the following features of the diagram.

The bands run between intersection points of singular lines for . These singular lines are, from (4.5), given by .

As the lattice sum is singular at while is not, the bands cannot cut singular lines except where they intersect.

As the band trajectory approaches an intersection point of singular lines, it follows that path along which the neighbouring singularities cancel.

There is a subtlety when singular lines are degenerate (the same trajectory of

*β*corresponds to two different pairs and ). This occurs for example at the lower right-hand side of figure 4, where the singularity lines and meet at . They both continue to the right as . It is clear from the band diagram that in fact the zero radius mode exactly follows the degenerate singular line.

In support of point (iv), we note that we can follow a mode with small *a* down to the singular line as . For example, if , and , the mode occurs at , above the singularity at . For , the mode has moved down to , while for , the mode is at . Thus, the mode tends down to the singular line, with its distance from the line scaling as . This has been indicated in figure 4 by drawing the mode line slightly displaced above the singularity line. Note also that the third band at the leftmost point in figure 4 joins up with the third band at the rightmost point.

In figure 5, we give the band diagram for the first three bands of a square array of circular holes with free edges, together with the singularity lines of , which we recall are just the dispersion curves for plane waves satisfying the Bloch–Floquet conditions. In contrast with the case of clamped boundaries, the diagram shows no band gaps within the frequency range studied. The behaviour of for the fundamental mode near is parabolic, as is typical for modes of plate structures (Graff 1975). Indeed, the band diagram for the perforated plate is similar to that for an unperforated plate, with bands tending to move down somewhat in frequency. Evidently, the band diagram of figure 5 is susceptible to a homogenization analysis, whereas that of figure 2 is not.

In figure 6, we compare the displacement patterns for modes calculated using the multipole approximation and the finite-element method, which agree well. Both methods show the expected displacement pattern of a perturbation of a plane wave associated with an array of circular holes. (Note that the phase variation of *W* across the unit cell is apparent in the plot of , not given here.)

## Acknowledgments

We acknowledge support from the EPSRC under grant EP/E035272/1 and from the Australian Research Council.

## Footnotes

- Received April 3, 2007.
- Accepted June 1, 2007.

- © 2007 The Royal Society