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 byEmbedded Imagewhere Embedded Image

Figure 1

An infinite plate with a doubly periodic array of circular voids and the elementary cell.

The out-of-plane elastic displacement Embedded Image satisfies the equation of motionEmbedded Image(2.1)where t is the time; h is the plate thickness; ρ is the mass density; Embedded Image 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 thatEmbedded Imagewhere (r, θ) are the polar coordinates and p is the angular frequency. In this case, equation (2.1) is transformed to the formEmbedded Image(2.2)where Embedded Image

The solution of (2.2) admits the representationEmbedded Image(2.3)whereEmbedded Image(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 representationEmbedded Image(2.5)with Wn defined byEmbedded Image(2.6)where, in (2.6), Jn is the Bessel function; Embedded Image is the Hankel function; and In and Kn are the modified Bessel functions related to Jn and Embedded Image, respectively, by (see Abramowitz & Stegun 1965, eqns 9.6.3 and 9.6.4)Embedded Image(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: Embedded Image. This means that the expansion covers the entire unit cell if Embedded Image. 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 Embedded Image and Embedded Image. 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 Embedded Image to that in the central cell Embedded Image asEmbedded Image(2.8)whereEmbedded Image(2.9)k0 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 Embedded Image, which are as follows.

  1. Clamped boundaryEmbedded Image(2.10)

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

Embedded Image(2.13)andEmbedded Image(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 k0. 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) implyEmbedded Image(3.1)andEmbedded Image(3.2)

Next, we address the case of the free boundary. The first condition in (2.11) together with (2.6) givesEmbedded Image(3.3)which (with the use of the standard Bessel equations for Embedded Image and Embedded Image, respectively) can be transformed to the formEmbedded Image(3.4)The second condition in (2.11) takes the formEmbedded Image(3.5)or equivalently,Embedded Image(3.6)Equation (3.6) can be simplified using standard Bessel equations, which giveEmbedded Image(3.7)andEmbedded Image(3.8)and similarly for the other two types of Bessel functions. Thus, we obtainEmbedded Image(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 Embedded Image. Their formula (9.1.79) isEmbedded Image(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 Embedded Image. We will also need the addition theorem for the modified Bessel functions K and I. This is obtained by replacing u,v,w by Embedded Image, and using the connection between Embedded Image and Embedded Image, and between Embedded Image and Embedded Image. The result isEmbedded Image(4.2)with a similar formula for In.

The Rayleigh identity for the function W(1), which satisfies the Helmholtz equation (2.4), is given byEmbedded Image(4.3)In (4.3), the quantity Embedded Image is the lattice sum of order l for the Hankel function H(1)Embedded Image(4.4)where Rq is defined in (2.9) and has the polar coordinate representation Embedded Image 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 Embedded Image (Chin et al. 1994), where various ways of evaluating the non-trivial lattice sums are reviewed by Moroz (2006). Here, we calculate Embedded Image using the representation obtained by repeated integration to accelerate convergence (Chin et al. 1994)Embedded Image(4.5)Here, Embedded Image denotes a pair of integers, which are summed over a sufficiently large range to guarantee the desired accuracy, while Embedded Image runs over vectors with the required quasi-periodicity of the Bloch vector k0. 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 Qh, 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 formEmbedded Image(4.6)whereEmbedded Image(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) isEmbedded Image(4.8)The version of this corresponding to acceleration by threefold integration isEmbedded Image(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:Embedded Image(4.10)andEmbedded Image(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 Embedded Image to Embedded Image

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

For the truncated system, arranging the unknown coefficients Embedded Image and Embedded Image into column vectors E and F, equations (5.1) and (5.2) can be combined into a single matrix equationEmbedded Image(5.3)where each of the four sub-blocks Embedded Image has the dimension Embedded Image

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

The equations for holes with clamped (5.3) or free (5.6) edges take the formEmbedded ImageThe dispersion equation then isEmbedded Imagewhich we solve numerically for a given value of k0 for the desired number of solutions β. The truncation order L we use here is 1, so the complex-valued matrix Embedded Image 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 Embedded Image: from Embedded Image to Embedded Image; from Embedded Image to Embedded Image; and (to the left) from Embedded Image to Embedded Image, with Embedded Image. The distance from the origin along the horizontal axis is the modulus of the Bloch vector k0. 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.

Figure 2

The band diagram giving Graphic as a function of Graphic for the case of clamped boundaries with hole radius Graphic.

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. 1995a,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 Embedded Image.

The determinant of the resulting 2×2 matrix isEmbedded Image(6.1)For the data of figure 2, with Embedded Image and Embedded Image, 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 Embedded Image) 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 Embedded Image and Embedded Image. 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 Embedded Image, 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.)

Figure 3

The fundamental mode displacement W as a function of position Graphic for clamped holes of radius Graphic, with Graphic: (a) analytical approximation and (b) numerical simulation in Comsol v. 3.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:Embedded Image(6.2)The band diagram corresponding to (6.2) is given in figure 4, together with lines marking the singularities of the lattice sum Embedded Image. Note the following features of the diagram.

Figure 4

The band diagram giving Graphic as a function of Graphic for the case of clamped holes in the zero radius limit. Solid lines, Bloch modes; dashed lines, singularities of the lattice sum Graphic.

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

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

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

  4. There is a subtlety when singular lines are degenerate (the same trajectory of β corresponds to two different pairs Embedded Image and Embedded Image). This occurs for example at the lower right-hand side of figure 4, where the singularity lines Embedded Image and Embedded Image meet at Embedded Image. They both continue to the right as Embedded Image. 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 Embedded Image. For example, if Embedded Image, Embedded Image and Embedded Image, the mode occurs at Embedded Image, above the singularity at Embedded Image. For Embedded Image, the mode has moved down to Embedded Image, while for Embedded Image, the mode is at Embedded Image. Thus, the mode tends down to the singular line, with its distance from the line scaling as Embedded Image. 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 Embedded Image, 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 Embedded Image for the fundamental mode near Embedded Image 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.

Figure 5

The band diagram giving Graphic as a function of Graphic for the case of holes of radius Graphic, with free edges, together with the dispersion curves for plane waves (dashed lines).

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 Embedded Image, not given here.)

Figure 6

The fundamental mode displacement W as a function of position Graphic for holes with free edges of radius Graphic, with Graphic: (a) analytical approximation and (b) numerical simulation in Comsol v. 3.2.


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


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


View Abstract