## Abstract

We consider wave propagation through a doubly periodic array of cavity cylinders in an isotropic elastic medium using the method of matched asymptotic expansions based on the assumptions that the scatterer size is much smaller than both the wavelength and the array periodicity. There is no restriction on the wavelength relative to the periodicity, and hence the method yields explicit approximations that describe the phenomena associated with periodic media. This is illustrated with square and hexagonal lattices.

## 1. Introduction

Wave propagation through periodic media has been a subject of extensive study for many years (see the online bibliography compiled by Dowling 2008). This paper is concerned with two-dimensional propagation in an elastic medium containing a doubly periodic lattice of circular–cylindrical inclusions. Accurate solution methods for this problem are provided by Poulton *et al.* (2000), Zalipaev *et al.* (2002) and Mei *et al.* (2003), and generalizations to oblique wave incidence and elliptical inclusions are described by Guenneau *et al.* (2003, 2005) and Guenneau & Movchan (2004). Low-frequency approximations for circular inclusions are obtained directly from their multi-pole formulation by Zalipaev *et al.* (2002), while other authors, for example Parnell & Abrahams (2006) and Andrianov *et al.* (2008), have used asymptotic homogenization to obtain low-frequency approximations to elastic wave propagation through periodic materials.

In the present work, the matched asymptotic expansions are used to obtain approximations for small scatterers that are perturbations of the quasi-periodic plane-wave solutions that exist in the absence of the scatterers. This method is used, for example, by Datta & El-Akily (1978) to solve the elastic wave-scattering problem. Here, the method follows that of McIver (2007), who investigated acoustic-wave propagation through a lattice of rigid scatterers. The main difference between the acoustic and elastic cases is that in the former only dilatational waves are present, while in the latter we need to take account of the coupling between shear and dilatational waves that arises from the boundary conditions. A notable feature of the modified matching procedure for the elastic case is that certain eigenfunctions must be included in the inner solution ahead of any obvious need for them. The main assumptions adopted are that, for waves with a typical wavenumber *k*, the scatterer size *a* is much smaller than both the scale *L* of the lattice periodicity and the wavelength 1/*k*. However, no restriction is placed on the wavelength relative to the periodicity *L* and hence it is possible to describe phenomena, such as stop bands, that are associated with the periodicity of the lattice. Craster *et al.* (2010) described a multiple scales approach to obtain results valid outside the low-frequency regime and near the edges of the Brillouin zone (although they do not apply their method to elastic waves). In contrast to the present work, application of the method of Craster *et al.* (2010) would require the numerical solution of a single-cell problem but, on the other hand, there would be no restriction on the size of an inclusion.

All solutions considered satisfy a Bloch condition that, for a specified value *β*_{0} of the Bloch wave vector ** β**, relates the solutions at corresponding points in different cells of the lattice. In the absence of the scatterers, and for the given Bloch wave vector

*β*_{0}, plane-wave solutions are possible for discrete values

*ω*

_{1},

*ω*

_{2},… of the frequency

*ω*. For each

*ω*

_{i}, there are

*M*≥1 plane waves corresponding to a particular pair (

*β*_{0},

*ω*

_{i})—these solutions may be shear or dilatational waves, or a mixture of the two. With the scatterers present, the asymptotic analysis yields an algebraic system to determine

*ω*for the

*M*perturbed modes that exist for each

**within a neighbourhood of (**

*β*

*β*_{0},

*ω*

_{i}) in the (

**,**

*β**ω*)-space. Explicit expressions for the frequencies are readily obtained that show how the mode frequencies depend on

**, the geometry, and the Lamé constants for the medium. Results are given to illustrate the appearance of local band gaps, the splitting and crossing of double modes, and the switching between dilatational and shear modes.**

*β*## 2. Problem formulation

We consider the propagation of waves through an isotropic elastic medium containing a doubly periodic lattice *Λ* of cavities, each of which is an infinitely long circular cylinder of radius *a*. Cartesian coordinates are chosen with origin *O* on the axis of one the cylinders, and with the *z*-axis directed along the axis of that cylinder; polar coordinates in the *x*–*y* plane with origin at *O* are denoted by (*r*,*θ*). In the *x*–*y* plane, scatterer *j* is associated with a local origin *O*_{j} located at the lattice point
2.1
for given independent vectors **a**_{1} and **a**_{2}.

The time-harmonic displacement of the medium outside the cavities is described by Navier's equation
2.2
where **u** is the displacement vector, *ω* is the angular frequency, *λ* and *μ* are Lamé constants and *ρ* is the density of the medium. Solutions of equation (2.2) are sought that satisfy the quasi-periodic ‘Bloch condition’
2.3
where ** β** is a prescribed wave vector, and traction-free boundary conditions are applied on the surface of each scatterer; thus, on

*r*=

*a*, 2.4 2.5 and 2.6 while equation (2.3) guarantees that the boundary conditions are applied throughout the lattice.

The main idea used here is to obtain perturbations of the plane-wave solutions that exist in the absence of the scatterers, and this is done using the method of matched asymptotic expansions. The wavenumbers for plane dilatational and shear waves are denoted, respectively, by *k*_{1} () and *k*_{2} () so that *k*_{2}=*Pk*_{1}, where *P*^{2}=(*λ*+2*μ*)/*μ*. It is assumed that the characteristic size *a* of a scatterer satisfies both *k*_{i}*a*≪1, *i*=1,2, and *a*/*L*≪1, where *L* is the length scale for the array periodicity, so that the scatterers are small relative to both the wavelengths and the array periodicity. The results given here are distinguished from those obtained through homogenization in that now each *k*_{i}*L* is allowed to be an order one quantity. This paper is concerned with waves that propagate in a direction normal to the axes of the cylinders, that is, in the *x*–*y* plane so that **u** does not depend on *z* (the same method is used by McIver (2007), to obtain results that apply to the simpler case of anti-plane shear waves).

By the Helmholtz representation,
2.7
where ** ψ**=(0,0,

*ψ*), so that 2.8 It follows that

*ϕ*and

*ψ*satisfy the two-dimensional Helmholtz equations 2.9 while the traction-free boundary conditions (2.5) and (2.6) on

*r*=

*a*become 2.10 and 2.11 In terms of

*ϕ*and

*ψ*, the Bloch conditions are 2.12 In particular, these conditions are satisfied by plane waves of the form 2.13 where

*β*_{m}=

**+**

*β***K**

_{m},

**=(**

*β**q*

_{1},

*q*

_{2})

^{T}is the prescribed Bloch vector, and each 2.14 is a reciprocal lattice vector with 2.15 The reciprocal lattice vectors have the property that, for any lattice vector

**R**

_{j}, 2.16 In this paper, results will be given for square and hexagonal lattices with reciprocal lattice vector pairs given respectively by

**b**

_{1}=

*L*

^{−1}(1,0),

**b**

_{2}=

*L*

^{−1}(0,1) and , , as illustrated in figure 1.

In the absence of the scatterers each of *ϕ*_{m} and *ψ*_{n} provides a solution to the Bloch problem provided *k*_{1}and *k*_{2} have values that ensure the field equations (2.9) are satisfied; in other words provided
2.17
where *β*_{m}=|*β*_{m}|. For example, for a square lattice of side *L* aligned with the coordinate axes, plane-wave solutions satisfying the Bloch conditions are
2.18
and
2.19
and the field equations are satisfied as long as
2.20
and
2.21
The results given here for the case when scatterers are present arise from the consideration of perturbations to combinations of such plane-wave solutions.

## 3. Solution by matched asymptotic expansions

The solutions obtained here are in terms of an outer solution, valid everywhere except for a neighbourhood of each scatterer, that matches with inner solutions valid only in the immediate vicinity of each scatterer. However, because of the Bloch conditions (2.12), the matching process used in the construction of solutions need only be carried out in a primary lattice cell chosen to be that containing the origin *O* of the global coordinates defined in the previous section. The particular scatterer associated with *O* is denoted by *C*. In addition to the global coordinates, local polar coordinates (*r*_{j},*θ*_{j}) are used that have origin at *O*_{j}.

To facilitate the solution, each lattice cell is divided into two overlapping regions. For the primary cell these are an outer region at distances *r*≫*a*, and an inner region within distances of the scatterer. A small parameter *ϵ*=*k*_{1}*a* is introduced, and in the inner region a scaled coordinate *ρ*=*r*/*a* is used. With these definitions, *k*_{1}*r*=*ϵρ* and *k*_{2}*r*=*Pϵρ*.

In the outer region, far from each scatterer, the solutions are constructed from solutions of the Helmholtz equations (2.9) that satisfy the Bloch conditions (2.12) and that are singular at the lattice points; such solutions are
3.1
and
3.2
where is associated with dilatational waves, and with shear waves. By Graf's addition theorem (Abramowitz & Stegun 1964, eqn 9.1.79)
3.3
and
3.4
where indicates summations over all , and the lattice sums are
3.5
Here *α*_{j} is the angle between the *x*-axis and **R**_{j} is measured in the anticlockwise direction, and the dashes indicate that **R**_{j}=0 is omitted from the summations. The lattice sums and in (3.5) have poles wherever, respectively, *k*_{1} and *k*_{2} have a value *β*_{m}, —see Linton (2010, eqn 3.18), for example—and these poles correspond to the plane-wave solutions discussed in §2. For a given ** β** and

*k*

_{i}(

*i*=1 or 2), the number of distinct vectors

*β*_{m}that have the same magnitude

*β*

_{m}is denoted by

*M*

_{i}≥1.

For square and hexagonal lattices, the locations of the poles of the lattice sums (or, equivalently, the plane waves that exist in the absence of the scatterers) are shown in figures 2 and 3 for values of the modulus *β* of the wave vector ** β** along the boundaries of the corresponding first irreducible Brillouin zones, as illustrated in figure 1. For all of the calculations in this paper, the Lamé constants are related by

*λ*/

*μ*=2.3 corresponding to a Poisson ratio

*ν*=0.35, and a non-dimensional frequency parameter

*ωL*/

*c*

_{2}≡

*k*

_{2}

*L*is used (the physical constants and frequency parameter are those used by Poulton

*et al.*(2000) and Zalipaev

*et al.*(2002)). In these figures it can be seen that, for some combinations of

**and the frequency**

*β**ω*, there are multiple plane-wave solutions and these can arise in two ways. First of all, as noted above, for one of the lattice sums there may be multiple distinct vectors

*β*_{m}with the same magnitude

*β*

_{m}; this can occur along lines as indicated by the two-pole curves in the figures, and also at isolated points. Second, multiple solutions can occur when the two lattice sums have poles corresponding to the same frequency, so that there are crossings of the curves for pure dilatational and pure shear waves. Results are given later for perturbations of all of these types of plane-wave solution, but full details of the derivation are given only for perturbations of purely dilatational waves (that is, in regions of the parameter space not close to the curves in figures 2 and 3 that correspond to shear waves).

Consider then the perturbation of dilatational plane waves with wavenumber *k*_{1}, so that the lattice sums have poles at *k*_{1}=±*β*_{m}, *m*=1..*M*_{1}, while the lattice sums for the shear waves are analytic functions of the frequency within neighbourhoods of these points. For each unique vector *β*_{m},
3.6
where is the area of one cell of the lattice, and the angles *τ*_{m} are defined through
3.7
see Linton (2010, equation 3.18). The lattice sums are written as follows:
3.8
and each is assumed to be an analytic function of *k*_{1} within neighbourhoods of each *k*_{1}=±*β*_{m}. Solutions are sought for *k*_{1} in a neighbourhood of *β*_{m} and it is assumed here that
3.9
where, for , *δ*_{m} is strictly of order one in *ϵ* as *ϵ*→0. This expression will be used within a neighbourhood of the points in (** β**,

*k*

_{1}) space that correspond to plane waves, so that the

*β*

_{m}and hence the

*δ*

_{m}may be distinct. It may be shown that the need to link the first appearance of singular terms in the outer solution of

*ϕ*with the non-singular leading-order outer solution for

*ϕ*requires

*f*(

*ϵ*)=

*ϵ*

^{2}(the matching would fail if this relation were incorrect) but, for simplicity, this will be adopted from the outset.

In view of equation (3.9), the matching may be carried out more conveniently, if the singular solutions of the Helmholtz equation defined in equation (3.1) are modified to be 3.10 where the plane-wave part 3.11 and the singular part 3.12

We will use *Φ* and to denote the outer and inner solutions for *ϕ*, and *Ψ* and to denote the outer and inner solutions for *ψ*. As the coupled boundary-value problem is homogeneous, the leading-order outer solution may be taken as strictly order zero in *ϵ* so that
3.13
where *Φ*^{(m)} denotes the outer solution up to terms in *ϵ*^{m}. (Because we consider only the perturbation of dilatational waves, shear waves do not appear in the leading-order outer solution.) In the following, *Φ*^{(m,l)} denotes the expansion up to *ϵ*^{l} of *Φ*^{(m)} after it is written in terms of the inner coordinate *ρ*. The inner solution up to terms in *ϵ*^{l} is denoted by , and denotes its expansion up to *ϵ*^{m} after it is written in terms of the outer coordinate *k*_{1}*r*. Matching is enforced by requiring for every *m* and *l* when both asymptotic forms are expressed in terms of the same coordinates (Crighton & Leppington 1973). The inner and outer expansions for the shear part of the solution are notated using the same principles.

In terms of the inner variable *ρ*=*r*/*a*, the field equations for the inner solutions and are
3.14
and
3.15
The inner solutions are constructed with the help of inner eigensolutions that each satisfy the Laplace equation, together with the homogeneous boundary conditions
3.16
and
3.17
It is straightforward to show that the required inner eigensolutions are given by the following pairs:
3.18

From equations (3.11) and (3.13), the inner expansion of the leading-order outer solution is
3.19
which indicates an initial inner development
3.20
where the term in *ν*_{11}(*ϵ*) is a possible intermediate term and, because and are coupled through the boundary conditions, there is a shear component to the inner problem so that
3.21
(higher order terms in the inner solution are dealt with later). Substituting (3.20) and (3.21) into the field equations (3.14) and (3.15) and equating the coefficients of the gauge functions in *ϵ*, we find that , , , , and are all harmonic functions that satisfy the homogeneous boundary conditions, and hence are constructed from the inner eigenfunctions detailed in equations (3.18).

The inner expansion (3.19) suggests that and might contain only constants. However, the appearance of terms in 2*θ* at order *ϵ*^{2} in the inner expansion (3.19) is significant. These terms could be matched with the inner solution simply by including appropriate eigenfunctions in at order *ϵ*^{2}, but it then proves impossible to satisfy the boundary conditions. The situation is resolved by including eigenfunctions in 2*θ* within and , in order to generate further terms at order *ϵ*^{2} in and through particular solutions of the field equations (3.14) and (3.15). In general, in order to satisfy the boundary conditions at order *ϵ*^{n+2}, the inner solution at order *ϵ*^{n} must contain singular eigenfunctions up to and including those in (*n*+2)*θ*.

With the above in mind, the leading-order inner solutions are written as
3.22
and
3.23
so that, in terms of the outer variables,
3.24
The outer expansion has terms no more singular than a quadrupole, and hence cannot be matched to any higher singularities associated in the leading-order outer solution; thus, for |*n*|>2. (It might be argued that more singular terms could be included in the outer solution, and then matched to the inner solution by including inner eigenfunctions at an appropriate order. However, this leads to unresolvable difficulties later.)

With possible intermediate terms included, the dilatational outer solution is continued as 3.25 so that, in particular, 3.26 The matching rule then gives 3.27

The right-hand equation in (3.24) shows that the leading term in the outer expansion of is at *ϵ*^{2}, and hence the leading-order outer shear solution is
3.28
(As we have assumed that the shear lattice sums are well behaved, there is no decomposition similar to equation (3.10) for the shear potentials. Also, because the leading inner term is strictly order one, when expressed in terms of the inner coordinates *Ψ*^{(2)} cannot be larger than order one as *ϵ*→0, and hence for |*n*|>2.) Thus, in particular,
3.29
and the matching rule gives
3.30
The inner expansion of the outer solution of *ϕ* up to order *ϵ*^{2} is
3.31
Thus, to match with the inner solution, *μ*_{11}(*ϵ*)=*ν*_{11}(*ϵ*), and the inner solutions must be continued as
3.32
and
3.33
(recall that at order *n* in the inner solution we must include singular eigenfunctions up to those in (*n*+2)*θ*). The potentials and satisfy the Poisson equations
3.34
and
3.35
as well as the boundary conditions (3.16) and (3.17). The solution forms needed to effect the matching are
3.36
and
3.37
where the ellipses indicate singular eigenfunctions in 3*θ* and 4*θ*. Further, the terms involving are particular solutions of (3.34) and (3.35), and are eigenfunctions, and the other terms are solutions of the Laplace equation to satisfy the boundary conditions and to effect the matching. From the boundary conditions
3.38
It follows that the outer expansions of and are, respectively,
3.39
and
3.40
The matching rules and yield the following:
3.41
3.42
3.43
and
3.44
To obtain a relationship between **u**_{1} and **u**_{I}, we need the order *ϵ*^{3} inner terms and which satisfy the Poisson equations
( and are given by terms in square brackets in, respectively, equations (3.32) and (3.33)). The required solutions of these Poisson equations are
and
where the ellipses indicate singular eigenfunctions in 4*θ* and 5*θ*. In these last equations, the terms in square brackets are particular solutions of the Poisson equations, and all other terms are eigenfunctions needed to satisfy the boundary conditions and effect the higher order matching. From the boundary conditions we get, in particular, **u**_{I}=−**u**_{1}/2.

Using equations (3.41) and (3.42), we now replace and in *B*_{0} (equation (3.27)), **u**_{0} (equation (3.43)) and **u**_{1} (equation (3.44)) to get
3.45
3.46
and
3.47
A system more amenable to further analysis is obtained by introducing
3.48
so that
3.49
and
3.50
and hence by substitution back into equation (3.48)
3.51
For a given ** β**, equation (3.51) provides an eigenvalue problem for the corresponding wave number

*k*

_{1}(which appears in each

*δ*

_{p}). The geometry of the lattice

*Λ*appears through the reciprocal lattice vectors in the definitions of each

*δ*

_{p},

*e*

_{1p}and

*e*

_{2p}.

For the perturbation of shear waves, a similar calculation yields
3.52
where , which is an eigenvalue problem for the shear wavenumber *k*_{2}. Comparing equations (3.51) and (3.52), we see that there is an extra term 1−*P*^{2} in (3.51). This arises from the form of the boundary condition for *σ*_{rr} in equation (2.10), which contains a non-derivative term in *ϕ*; this generates a monopole, and hence the additional term, in the perturbed dilatational wave, but not in the perturbed shear wave.

It is also possible for the frequency of an unperturbed dilatational wave to coincide with the frequency of an unperturbed shear wave. The simultaneous perturbation of dilatational and shear waves will, in general, result in a pair of coupled eigenvalue problems for the wavenumbers *k*_{1} and *k*_{2}, namely
3.53
and
3.54
where *p*_{1}=1..*M*_{1}, *p*_{2}=1..*M*_{2}. Here the terms with the bar refer to the perturbed dilatational wave, and the terms with a hat refer to the perturbed shear wave. It is noteworthy that in the case *M*_{1}=*M*_{2}=1, that is, the perturbation of one dilatational and one shear wave, the off-diagonal terms in the system matrix are zero and the equations uncouple. Equations (3.53) and (3.54) include the cases of the perturbed dilatational waves and shear wave already considered separately; with , equation (3.53) reduces to (3.54), and with equation (3.54) reduces to (3.52).

## 4. Results

In this section, we give some examples of explicit approximations to the dispersion relation for perturbed dilatational and shear waves, obtained from the eigenvalue problems in equations (3.51)–(3.54) with the aid of the computer algebra package Mathematica. In addition, some of these approximations are compared with numerical results for *λ*/*μ*=2.3, so that *P*^{2}=(*λ*+2*μ*)/*μ*=4.3. Results are given on the boundary of the irreducible region of the first Brillouin zone (labelled GMK in each part of figure 1) for square and hexagonal lattices. The areas of the corresponding primitive cells are, respectively, and .

### (a) Perturbation of one plane wave

For *M*=1, equation (3.51) for the perturbation of a dilatational wave is
4.1
while equation (3.52) for the perturbation of a shear wave gives
4.2
These expressions give perturbations of any of the one-pole solutions illustrated in the dispersion diagrams in figures 2 and 3, regardless of their frequency. As noted at the end of the previous section, this includes any points in the diagrams where one-pole dilatational- and shear-wave modes cross. For the lowest dilatational mode and two lowest shear modes, and along the line GM of the irreducible Brillouin zone, the present approximations are compared in figure 4 with numerical calculations made using the multi-pole method described by Zalipaev *et al.* (2002), and generally a good agreement is observed. The divergence of the curves near the ends of GM is because of the presence of higher order poles—see figures 2 and 3. Zalipaev *et al.* (2002) give asymptotic results that differ from equations (4.1) and (4.2) because, in their calculations, the lattice sums were taken as order one when, as can be seen from equations (3.8) and (3.9), they are singular in the limit *ϵ*→0; further details are given by Guo (2011).

### (b) Local band gaps

Approximations at appropriate isolated points in the dispersion diagram lead to estimates of local band gaps, and this is illustrated here by consideration of a two-pole solution corresponding to the perturbation of shear waves at the point M in the reciprocal lattice. For the square lattice, the unperturbed values considered are (*q*_{1}*L*,*q*_{2}*L*,*k*_{2}*L*)=(*π*,0,3*π*), and within some neighbourhood of this point the appropriate forms for the *β*_{m} are
4.3
(this is because, when (*q*_{1}*L*,*q*_{2}*L*)=(*π*,0), *β*_{1}*L*=*β*_{2}*L*=3*π*≡*k*_{2}*L*). For *q*_{2}=0 and in the limit *q*_{1}*L*→*π*, the positive roots of equation (3.52) are
4.4
Similarly, for the hexagonal lattice, with there are again two poles and the appropriate forms for the *β*_{m} are
4.5
For *q*_{1}*L*=0, and as , the positive roots of equation (3.52) are
4.6
Equations (4.4) and (4.6) illustrate explicitly the appearance of local band gaps as the radius of the cylinder is increased from zero. In each case, the upper point of the local band gap is independent of the Lamé constants and is determined by the geometry alone.

### (c) Mode splitting and mode switching

Two-pole lines in the unperturbed dispersion diagrams in figures 2 and 3 will, in general, split into two modes when the scatterers are present. This is shown first for the lowest dilatational two-pole mode on MK with the hexagonal lattice, for which and *q*_{1}*L*∈(0,2*π*/3), and the appropriate forms for *β*_{m} are
4.7
The corresponding positive roots of (3.51) are
4.8
and
4.9
where
4.10
and
4.11
These solutions are graphed in figure 5, which shows how the gap increases with the radius of the cylinder. In figure 3, it may be seen that a two-pole shear mode crosses near K but, for simplicity, this is not taken into account here.

Split modes may also cross, and this is illustrated using the lowest (shear) two-pole line on MK in figure 2 for the square lattice, for which *q*_{1}*L*=*π* and *q*_{2}*L*∈(0,*π*), and the appropriate forms for *β*_{m} are
4.12
For these values, the positive roots of equation (3.52) are
4.13
and
4.14
where
The gap between the two curves depends on the radius of the scatterer and is
4.15
so that, to a first approximation, the two lines cross at
4.16
to the left of this point (*k*_{2}*L*)_{1} is the upper mode, and (*k*_{2}*L*)_{2} is the lower mode. The above approximations are compared with numerical solutions in figure 6*a*. The approximations are worst near the right-hand side of the range of *q*_{2}*L* owing to the proximity of the four-pole point at K. One feature of the present method is that solutions valid in the neighbourhood of higher order poles blend smoothly into lower order solutions, and hence the former may be used outside their apparent range of validity. This may be shown explicitly in simple cases, and is illustrated graphically in figure 6*b*, which uses the four-pole approximation (not given explicitly here) from the point K across the whole of MK.

As well as mode splitting, there can also be mode switching where, in the unperturbed dispersion diagram, shear and dilatational modes cross. This is illustrated in figure 7 for the intersection of the lowest two-pole shear modes in KG with a dilatational mode. The local solution at the three-pole intersection point correctly shows the switching between the dilatational mode and one of the shear modes. The approximation degrades towards the ends of the range because of the proximity of higher order poles.

## Acknowledgements

We are grateful to Dr Ian Thompson for allowing us to use his Fortran code for the calculation of lattice sums.

- Received January 25, 2011.
- Accepted April 21, 2011.

- This journal is © 2011 The Royal Society