\begin{opening} \title{The Linear Wave Response of a Periodic Array of Floating \\ Elastic Plates}

C. D. \surname{Wang}\email{cvewcd@nus.edu.sg} \institute{Department of Civil Engineering, National University of Singapore, Kent Ridge, Singapore 119260, Singapore.} M. H. \surname{Meylan}\email{meylan@math.auckland.ac.nz} \institute{Department of Mathematics, University of Auckland, Private Bag 92019, Auckland, New Zealand.} R. \surname{Porter}\email{Richard.Porter@bristol.ac.uk} \institute{School of Mathematics, University of Bristol, Bristol, BS8 1TW, UK.} \institute{}

\runningtitle{Response of a Periodic Array of Plates} \runningauthor{Wang, Meylan \& Porter}

The problem of an infinite periodic array of identical floating elastic plates subject to forcing from plane incident waves is considered. This study is motivated by the problem of trying to model wave propagation in the marginal ice zone, a region of ocean consisting of an arbitrary packing of floating ice sheets. It is shown that the problem considered can be formulated exactly in terms of the solution to an integral equation in a manner similar to that used for the problem of wave scattering by a single elastic floating plate, the key difference here being the use of a modified periodic Green function. The convergence of this Green function in its original form is poor, but can be accelerated by a transformation. It is shown that the results from the method satisfy energy conservation and that in the particular case of a fixed rigid rectangular plate which spans the periodicity uniformly the solution reduces to that for a two-dimensional rigid dock. We also present solutions for a range of elastic plate geometries.

\keywords{Ice sheet, elastic plates, periodic Green function, water waves}

\end{opening}

Introduction

The thin elastic plate of shallow draft floating on the surface of water can be used to model a range of physical systems, for example ice floes or so-called very large floating structures, which are of particular interest in proposals to build offshore floating runways. It is also of theoretical importance as one of the simplest models of hydroelasticity. It is not surprising therefore that the floating elastic plate has been the subject of a significant body of research. Much of this research, especially that which was motivated by the construction of very large floating structures, is summarised in the review papers by Kashiwagi Kashiwagi00 and Watanabe, Utsunomia, and Wang Watanabe04. The review paper by Squire et al. Squire_Review summarises the research prior to 1995 which was motivated by modelling sea ice floes, but does not include the more recent work, Meylan and Squire jgrfloecirc and Meylan JGR02, in which the solution for wave scattering by a three-dimensional ice floe is presented.

Our aim here is to present a method of solution to the problem of wave scattering by an infinite periodic array of floating elastic plates of finite dimension. This study is motivated by trying to further understanding of the problem of wave scattering in the marginal ice zone. The marginal ice zone consists of a vast field of broken ice floes each of which have a thickness (typically [math]\displaystyle{ \approx 1 }[/math]-[math]\displaystyle{ 2 }[/math]m) much smaller than the typical horizontal lengthscale of the floe, or of the wavelength of waves supported by the floe, and it is usual to model the ice floe as a thin elastic plate. These ice floes are subject to intense wave forces whose source originates from the open ocean. On account of the fact that the floes are themselves able to support wave propagation, wave energy is capable of travelling large distances within the marginal ice zone where it can assist ice breakup.

Previous work aimed at trying to model wave interaction by isolated three-dimensional regular and irregular ice floes has been performed in Meylan and Squire jgrfloecirc and Meylan JGR02, respectively. For more complicated arrangements involving an arbitrary, but finite, number of irregular ice floes see Peter and Meylan JFM_malte04, although the number of floes that can be dealt with is limited by numerical considerations. Although this latter work includes the effect of interaction between neighbouring ice floes, it does not allow one to capture the effects of a plane wave incident upon a very large linear array of ice floes. One method of overcoming this, which also allows significant analytic progress to be made is to consider wave scattering by an infinite periodic array of identical ice floes, being irregular in shape. Of course ice fields are not periodic, but this arrangement does provide some information especially if some kind of statistical averaging of the results over the floe geometries and spacings is performed. For example, we might expect that a random array of ice floes would scatter equivalently to an average of periodic arrays of identical floes of different spacings and geometries. Furthermore, the solution readily lends itself to an assumption that the rows of periodic ice floes are wide-spaced. This method was used by Evans and Porter Evans95 and would allow accurate approximations to wave scattering by multiple rows of non-identical ice floes. Central to such a procedure is knowledge of the scattering process for each row in isolation.

The problem of determining the scattering of waves by periodic arrays of obstacles subject to wave forcing has received considerable research attention and spans a broad range of physical disciplines including solid-state physics, acoustics, optics, etc. In many applications, the interest centres on arrangements which are periodic in two directions (for example, the study of crystallography).

In the present context of water wave propagation and its interaction with flexible surface structures, Chou Chou98 has investigated the effect of an infinite doubly-periodic array of elastic plates on wave propagation. In the modelling of the plate equations, Chou incorporates both bending stiffness and tension effects, so that the discussion of the results not only includes the case of pure bending of elastic plates in the absence of compression forces (as considered here), but also, by setting the stiffness to zero, pure tensional effects which would describe, for example, periodic arrays of taught membranes. However, in contrast to the work described here, the doubly-periodic configuration allows significant simplification in the solution procedure by applying Floquet's theorem to reduce the problem to one on a finite domain with periodic boundary conditions. Moreover, problems involving infinite doubly-periodic structures only offer information about the possibility of wave propagation throughout the array (in the form of so-called pass-bands or stop-bands) and cannot address the diffraction of plane waves from infinity.

For arrays which are periodic in one direction only, the situation is different and diffraction grating effects occur. Thus, for an incident plane wave of a particular given wave frequency, a finite number of distinct plane waves propagating away from the array at certain discrete angles will occur. In the context of water waves and fixed periodic arrays, Twersky Twersky52 was able to solve the problem of a periodic array of vertical circular cylinders. The uniformity of the configuration in the depth coordinate implies that the resulting equations also describe two-dimensional acoustic wave scattering, in this case by circular cylinders. The problem of Twersky Twersky52, who used Schl\"{o}milch series to sum slowly-convergent series involving Hankel functions, was re-considered by Linton and Evans Linton_evans93 who used a so-called multipole method. For periodic arrays of rectangular cylinders extending uniformly through the depth, Fernyhough and Evans Fernyhough95, used domain decomposition and mode matching to derive an integral equation formulation to the problem. In order to consider more general cylinder profiles, boundary integral methods are inevitable and require the use of a periodic Green function. In its most basic form, the Green function consists of a series involving Hankel functions which is slowly convergent and unsuitable for numerical computation. Hence Linton Linton98 compared a number of different representations of the periodic Green function, designed to increase the convergence characteristics. Porter and Evans Porter_evans99 used the work of Linton Linton98 to compute so-called Rayleigh-Bloch waves (or trapped waves) along a periodic array of cylinders of arbitrary cross-section. A number of papers in recent years have concentrated on similar ideas, to those of Porter and Evans Porter_evans99, a primary motivation being the connection between large wave responses in large finite arrays of cylinders and the trapped waves in infinite periodic arrays (see Maniar and Newman Maniar_newman97).

In contrast to the body of work cited in the previous paragraph, the present work differs significantly in two respects. First, the problem involving periodic elastic plates on the water surface implies that, though the periodicity extends only in one direction, the problem uses fully three-dimensional plates. Secondly, the scattering obstacles are neither fixed (as above) nor rigid, but are in fact themselves capable of supporting wave motions. Thus the coupling between the plate motion and the motion of the water adds an extra degree of complication to the problem. We consider a special Green function which is build up out of multiple images as was done in Kashiwagi kashiwagi91 for the case of wave scattering in a wave tank with rigid walls. Our method may be considered as generalisation of Kashiwagi in which we allow the waves to be incident at an angle and the body spacing to be arbitrary. We also use a different method to sum the infinite series of Green functions that was used by Kashiwagi.

The outline of the paper is as follows. In section 2 the non-dimensional small-amplitude equations for time-harmonic motion are formulated. In section 3, an integral equation is derived connecting the motion of the plate to that of the water. In doing so the periodic Green function is defined in terms of the three-dimensional free-surface Green function. In section 4, a variational principle for the motion of the plate is defined and a Rayleigh-Ritz approximation is used to reduce it to a linear system of equations. The same approximation is used in the integral equation to provide a second set of linear equations coupling the motion of the plate to the water. Section 5 concentrates on improving the otherwise slow convergence of the periodic Green function and section 6 describes how the far-field reflected and transmitted waves are calculated from the integral formulation. Finally, in section 7 a range of numerical results are presented and discussed and the work is concluded in section 8.

Problem Formulation: An Infinite Array of Elastic Plates

We begin by formulating the problem. Cartesian coordinates [math]\displaystyle{ (x,y,z) }[/math] are chosen with [math]\displaystyle{ z }[/math] vertically upwards such that [math]\displaystyle{ z=0 }[/math] coincides with the mean free surface of the water. An infinite array of identical thin elastic plates float on the surface of the water, periodically spaced along the [math]\displaystyle{ y }[/math]-axis with uniform separation [math]\displaystyle{ l }[/math]. The problem is to determine the motion of the water and the plates when plane waves are obliquely-incident from [math]\displaystyle{ x=-\infty }[/math] upon the periodic array of plates.

The plates are assumed to be of zero draft and occupy [math]\displaystyle{ {\bf x} \in \Delta_m\lt math\gt , }[/math]-\infty < m < \infty[math]\displaystyle{ on }[/math]z=0[math]\displaystyle{ where }[/math]{\bf x} = (x,y)</math> represents the Cartesian vector lying in the mean free surface. The plates are assumed to have arbitrary shape, and periodicity implies that if [math]\displaystyle{ {\bf x} \in \Delta_0 }[/math], then [math]\displaystyle{ {\bf x}_m = (x,y+ml) \in \Delta_m }[/math], [math]\displaystyle{ -\infty \lt m \lt \infty }[/math]. This array of plates is shown in Figure (fig_array).

Elastic plate equations

The equation of motion for the plate elevation [math]\displaystyle{ W(x,y,t) }[/math], where [math]\displaystyle{ t }[/math] is time, is given by the thin elastic plate (or Kirchhoff) equation,

[math]\displaystyle{ D \nabla_h^{4} W + \rho_{i} h \frac{\partial^{2}W}{\partial t^{2}} = p, (floemot) }[/math]

\cite[pp. 79--104]{Timoshenko_woinowsky-krieger59} where [math]\displaystyle{ D }[/math] is the flexural rigidity, [math]\displaystyle{ \rho_i }[/math] is the density, and [math]\displaystyle{ h }[/math] the uniform thickness of the elastic plate. Here, [math]\displaystyle{ p }[/math] is the excess pressure (i.e. excluding atmospheric pressure and the weight of the plate) exerted by the fluid and [math]\displaystyle{ \nabla_h^2 }[/math] is the two-dimensional horizontal Laplacian in the plane [math]\displaystyle{ z=0 }[/math] (we shall reserve [math]\displaystyle{ \nabla^2 }[/math] to mean the three-dimensional Laplacian). In addition, on the free edges of the plates, [math]\displaystyle{ {\bf x} \in \partial \Delta_m }[/math], boundary conditions expressing the vanishing of bending moment and shearing stress apply, which are written as

[math]\displaystyle{ (boundary1) \left[ \nabla_h^2 - (1-\nu) \left(\frac{\partial^2}{\partial s^2} + \kappa(s) \frac{\partial}{\partial n} \right) \right] W = 0, }[/math] [math]\displaystyle{ (boundary2) \left[ \frac{\partial}{\partial n} \nabla_h^2 +(1-\nu) \frac{\partial}{\partial s} \left( \frac{\partial}{\partial n} \frac{\partial}{\partial s} -\kappa(s) \frac{\partial}{\partial s} \right) \right] W = 0, }[/math]

where [math]\displaystyle{ \nu }[/math] is Poisson's ratio and

[math]\displaystyle{ \nabla_h^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} = \frac{\partial^2}{\partial n^2} + \frac{\partial^2}{\partial s^2} + \kappa(s) \frac{\partial}{\partial n}. }[/math]

Here, [math]\displaystyle{ \kappa(s) }[/math] is the curvature of the boundary, [math]\displaystyle{ \partial \Delta_m }[/math], as a function of arclength [math]\displaystyle{ s }[/math] along [math]\displaystyle{ \partial \Delta_m }[/math]; [math]\displaystyle{ \partial/\partial s }[/math] and [math]\displaystyle{ \partial/\partial n }[/math] represent derivatives tangential and normal to the boundary [math]\displaystyle{ \partial \Delta_m }[/math].

Equations ((boundary1)) and ((boundary2)) can be established from Porter and Porter Porter_porter04 who considered the more complicated case of a plate of varying thickness and derived the edge conditions from a variational principle similar to the one used later in this paper. A direct derivation of the edge conditions from the underlying constitutive equations for a thin elastic plate can be found elsewhere, for example \cite[Chapter 8]{Veubeke79}. We remark that in the case where the boundary of the elastic plates is piecewise linear, such that [math]\displaystyle{ \kappa(s)=0 }[/math], the equations above reduce to the simpler form

[math]\displaystyle{ \left( \frac{\partial^2}{\partial n^2} + \nu \frac{\partial^2}{\partial s^2} \right) W = 0 }[/math]

and

[math]\displaystyle{ \left( \frac{\partial^3}{\partial n^3} + (2-\nu) \frac{\partial^3}{\partial n \partial s^2} \right) W = 0. }[/math]

Equations of motion for the water

We consider water of infinite depth. Under the usual assumptions of small-amplitude waves above an incompressible and irrotational fluid, there exists a velocity potential [math]\displaystyle{ \Phi(x,y,z,t) }[/math] satisfying Laplace's equation together with the appropriate linearised boundary conditions, namely

[math]\displaystyle{ \left. \begin{matrix}[c]{rcl} \nabla^2 \Phi & = & 0,\qquad -\infty\lt z\lt 0, \\ \noalign{\vskip4pt} \Phi, \; |\nabla \Phi| & \to & 0, \qquad z \to -\infty, \\ \noalign{\vskip4pt} {\displaystyle \frac{\partial \Phi}{\partial z}} & = & {\displaystyle \frac{\partial W}{\partial t}},\qquad z=0, \\ \noalign{\vskip4pt} {\displaystyle -\rho \frac{\partial\Phi}{\partial t} - \rho gW} & = & p,\qquad z=0, \end{matrix} \right\} (watermot) }[/math]

where the pressure, [math]\displaystyle{ p }[/math], has already been introduced in ((floemot)), [math]\displaystyle{ \rho }[/math] is the water density, and [math]\displaystyle{ g }[/math] is gravitational acceleration. [math]\displaystyle{ W }[/math] is the time-dependent displacement of the plate when [math]\displaystyle{ \mathbf{x} \in \Delta_m\lt math\gt . When }[/math]\mathbf{x} \not\in \Delta_m</math>, we take [math]\displaystyle{ p=0 }[/math] in the above and now [math]\displaystyle{ W }[/math] represents the elevation of the water surface. In addition to ((watermot)), radiation conditions at infinity need to be applied and these will formally be introduced later.

Non-dimensionalising the variables

We non-dimensionalise the spatial variables with respect to a length parameter [math]\displaystyle{ L }[/math] (for example, [math]\displaystyle{ L }[/math] may be derived from the area of the plate or [math]\displaystyle{ L }[/math] may be the characteristic length [math]\displaystyle{ \left( D/\rho g \right)^{1/4} }[/math]) and the time variables with respect to [math]\displaystyle{ \sqrt{L/g.} }[/math] The dimensionless variables are therefore given by

[math]\displaystyle{ (\bar{x},\bar{y},\bar{z},\bar {W}) = \frac{1}{L}(x,y,z,W),~ \bar{t} = t\sqrt{\frac{g}{L}},~ \bar{p} = \frac{p}{\rho g L} \quad \mbox{and} \quad \bar{\Phi} = \frac{\Phi}{L\sqrt{L\,g}}. }[/math]

Using the dimensionless variables equation ((floemot)) combined with the last equation of ((watermot)) becomes

[math]\displaystyle{ \beta \bar{\nabla}_h^{4} \bar{W} + \gamma \frac{\partial^{2}\bar{W}}{\partial\bar{t}^{2}} = \bar{p} = -\left. \frac{\partial\bar{\Phi}}{\partial\bar{t}} \right|_{z=0} -\bar{W}, (dimless_floepressmot_1) }[/math]

where the dimensionless parameters associated with the motion of the plate, [math]\displaystyle{ \beta }[/math] and [math]\displaystyle{ \gamma }[/math], representing the `stiffness' and `mass loading' of the plate respectively, are given by

[math]\displaystyle{ \beta = \frac{D}{\rho g L^{4}} \quad \mbox{and} \quad \gamma = \frac{\rho_{i}h}{\rho L}. }[/math]

This notation is based on Tayler \cite[pp. 122--128]{Tayler86}.

Hereafter, we will work with dimensionless variables only but omit the overbar from all variables for reasons of clarity. Note that from now on we will use [math]\displaystyle{ l }[/math] to mean the non-dimensional floe separation [math]\displaystyle{ \bar{l} = l/L }[/math].

The single frequency equations

We will consider the solution for a single frequency and we can therefore represent the displacement and the potential as the real parts of complex functions in which the time dependence is [math]\displaystyle{ \exp^{-i\omega t} }[/math] where [math]\displaystyle{ \omega }[/math] is the dimensionless radian frequency, i.e.

[math]\displaystyle{ \begin{matrix} W( x,y,t) & = & \mbox{Re} \left[ \left(\frac{i}{\omega}\right) w(x,y) \exp^{-i\omega t} \right],\\ \Phi(x,y,z,t) & = & \mbox{Re} \left[ \phi(x,y,z) \exp^{-i\omega t}\right], \end{matrix} }[/math]

where we have introduced an additional scaling for [math]\displaystyle{ W }[/math] to simplify the equations which follow.

Therefore equation ((dimless_floepressmot_1)) becomes

[math]\displaystyle{ \beta \nabla_h^{4} w(x,y) + (1 - \omega^{2} \gamma) w(x,y) = \omega^2 \phi(x,y,0). (maineqn) }[/math]

Combining the non-dimensionalisation with the time assumption we also have from ((watermot))

[math]\displaystyle{ \frac{\partial\phi}{\partial z} = w, \qquad z=0 (kineqn) }[/math]

with

[math]\displaystyle{ \left. \begin{matrix}[c]{rcl} \nabla^{2}\phi & = & 0,\qquad -\infty \lt z \lt 0,\\ \noalign{\vskip4pt} \phi,\; |\nabla \phi| & \to & 0, \qquad z \to -\infty,\\ \noalign{\vskip4pt} \end{matrix} \right\} (watermot2) }[/math]

For [math]\displaystyle{ \mathbf{x} \not\in \Delta_m }[/math], ((maineqn)) still holds, but with [math]\displaystyle{ \beta = \gamma = 0 }[/math] and combining with ((kineqn)) gives the usual free-surface condition

[math]\displaystyle{ \frac{\partial\phi}{\partial z} - k \phi = 0, \qquad \mbox{on \lt math\gt z=0 }[/math]}, \qquad \mbox{where [math]\displaystyle{ k = \omega^2 }[/math]} </math>

and [math]\displaystyle{ k }[/math] is the dimensionless wavenumber (i.e. the dimensionless wavelength is [math]\displaystyle{ \lambda=2\pi/k }[/math]). The velocity potential also satisfies the Sommerfeld radiation condition as [math]\displaystyle{ |\mathbf{x}| \rightarrow\infty }[/math],

[math]\displaystyle{ \lim_{|\mathbf{x}| \rightarrow \infty} \sqrt{|\mathbf{x}| }\left( \frac{\partial}{\partial |\mathbf{x}| } - i k \right) \left( \phi-\phi^{\rm in} \right) = 0, (sommerfeld) }[/math]

(see Wehausen and Laitone Weh_Lait). In equation ((sommerfeld)), [math]\displaystyle{ \phi^{\rm in} }[/math] is the incident wave potential given by

[math]\displaystyle{ \phi^{{\rm in}} = \frac{A}{k} \exp^{ik (x\cos\theta+y\sin\theta)}\,\exp^{kz}, (phi_inc) }[/math]

where [math]\displaystyle{ A }[/math] is the dimensionless amplitude and [math]\displaystyle{ \theta }[/math] is the direction of propagation of the wave.

Transformation to an Integral Equation

We now Floquet's theorem (Scott Scott98) (also called {\it the assumption of periodicity} in the water wave context) which states the displacement from adjacent plates differ only by a phase factor. If the potential under the {\em central} plate [math]\displaystyle{ \Delta_{0} }[/math] is given by [math]\displaystyle{ \phi( \mathbf{x}_{0},0) }[/math], [math]\displaystyle{ \mathbf{x}_{0}\in\Delta_{0} }[/math], then by Floquet's theorem the potential satisfies

[math]\displaystyle{ \phi(\mathbf{x}_{m},0) = \phi(\mathbf{x}_{0},0) \exp^{im\sigma l}, (floq1) }[/math]

and the displacement of the plate [math]\displaystyle{ \Delta_{m} }[/math] satisfies

[math]\displaystyle{ w(\mathbf{x}_{m}) = w(\mathbf{x}_{0}) \exp^{im\sigma l}, (floq2) }[/math]

where [math]\displaystyle{ \mathbf{x}_{m} \in \Delta_{m} }[/math], [math]\displaystyle{ -\infty \lt m \lt \infty }[/math] and the phase difference is [math]\displaystyle{ \sigma = k\sin\theta }[/math] (see, for example, Linton Linton98).

A standard approach to the solution of the equations of motion for the water ((maineqn)), ((kineqn)), ((watermot2)) is to transform these equations into a boundary integral equation using the free-surface Green function for infinite depth (see Weh_Lait,Kim65). In doing so we obtain

[math]\displaystyle{ \phi(\mathbf{x},0) = \phi^{\rm in} (\mathbf{x},0) +\sum_{m=-\infty}^{\infty} \int_{\Delta_{m}} G(\mathbf{x},0 ;\mbox{\boldmath\lt math\gt \xi }[/math]}) \left[ k\phi(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]},0)

- w (\mbox{\boldmath[math]\displaystyle{ \xi }[/math]})  \right]  d\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}
 (phi_w_m)

</math>

where [math]\displaystyle{ \mbox{\boldmath }[/math]\xi[math]\displaystyle{ } = (\xi,\eta) }[/math] and [math]\displaystyle{ G(\mathbf{x},z;\mbox{\boldmath }[/math]\xi[math]\displaystyle{ }) }[/math] is the free-surface Green function satisfying

[math]\displaystyle{ \left. \begin{matrix}[c]{rcl} \nabla^2 G & = & 0,\qquad -\infty \lt z \lt 0, \\ \noalign{\vskip4pt} \displaystyle{\frac{\partial G}{\partial z}} - k G & = & -\delta(\mathbf{x}-\mbox{\boldmath\lt math\gt \xi }[/math]}),\qquad z=0,

\\ \noalign{\vskip4pt}
G, \; |\nabla G| & \to & 0,\qquad z \to -\infty.

\end{matrix} \right\}

</math>

which, on [math]\displaystyle{ z=0 }[/math], is given by

[math]\displaystyle{ \begin{matrix} G(\mathbf{x},0;\mbox{\boldmath\lt math\gt \xi }[/math]}) & = &

-\frac{1}{4\pi} \Bigg( \frac {2}{\left|

\mathbf{x}-\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}\right|}

-\pi k \Big[ \mathbf{H}_{0}\left( k
\left| \mathbf{x}-\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}\right| \right)
 \\
& & \quad \quad \quad
+Y_{0}\left( k\left| \mathbf{x}-\mbox{\boldmath[math]\displaystyle{ \xi }[/math]} \right|
\right) -2i\pi J_{0}\left( k
\left| \mathbf{x}-\mbox{\boldmath[math]\displaystyle{ \xi }[/math]} \right| \right)  \Big] \Bigg),
 (green)

\end{matrix}</math>

In the above [math]\displaystyle{ \mathbf{H}_{0} }[/math] is the Struve function of order zero, [math]\displaystyle{ J_{0} }[/math] is the Bessel function of the first kind of order zero and [math]\displaystyle{ Y_{0} }[/math] is the Bessel function of the second kind of order zero (Abramowitz and Stegun \cite[Chapters 9 and 12]{abr_ste}).

Using ((floq1)) and ((floq2)) in ((phi_w_m)) we obtain

[math]\displaystyle{ \begin{matrix} (bem_eq_1) \phi (\mathbf{x},0) & = & \phi^{\rm in}(\mathbf{x},0) \\ & & \! +\!\! \sum_{m=-\infty}^{\infty} \int_{\Delta_{0}} \!\! G (\mathbf{x},0 ;\mbox{\boldmath\lt math\gt \xi }[/math]}+(0,ml)) \left[ k\phi(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]},0)

- w(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}) \right] \exp^{im\sigma l}
\,d\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}

\end{matrix}</math>

which can be written alternatively as

[math]\displaystyle{ (bem_eq_2) \phi(\mathbf{x},0) = \phi^{\rm in}(\mathbf{x},0) +\int_{\Delta_{0}} G_{\mathbf{P}} (\mathbf{x} ;\mbox{\boldmath\lt math\gt \xi }[/math]}) \left[ k\phi(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]},0)

- w(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]})  \right] \,d\mbox{\boldmath[math]\displaystyle{ \xi }[/math]},

</math>

where the kernel [math]\displaystyle{ G_{\mathbf{P}} }[/math] (referred to as the {\it periodic Green function}) is given by

[math]\displaystyle{ G_{\mathbf{P}} (\mathbf{x};\mbox{\boldmath\lt math\gt \xi }[/math]})

= \sum_{m=-\infty}^{\infty} G \left(\mathbf{x},0;\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}+(0,ml)
\right) \exp^{im\sigma l}.
 (G_periodic)

</math>

Solution of the Integral Equation

The integral equation ((bem_eq_2)) is identical to the integral equation for a single plate (see Meylan JGR02) except for the modification to the Green function. Furthermore, the periodic Green function [math]\displaystyle{ G_{\mathbf{P}} }[/math] has the same singularity at [math]\displaystyle{ \mathbf{x}= \mbox{\boldmath\lt math\gt \xi }[/math]}</math> as the standard free-surface Green function [math]\displaystyle{ G }[/math]. We solve equation ((bem_eq_2)) using a higher-order method, which is explained in detail in Wang and Meylan Wang04. We will briefly outline the solution method here.

We expand the plate potential and displacement as

[math]\displaystyle{ \left. \begin{matrix}{c} {\displaystyle w(\mathbf{x})\approx \sum_{i=1}^{3q}w_{i}\chi _{i}(\mathbf{x}) = \vec{\chi}^{T}(\mathbf{x})\vec{w} } \\ {\displaystyle \phi (\mathbf{x},0)\approx \sum_{i=1}^{3q}\phi_{i}\chi_{i}(\mathbf{x}) =\vec{\chi}^{T}(\mathbf{x})\vec{\phi} } \end{matrix} \right\} (expansion) }[/math]

where [math]\displaystyle{ \vec{w} }[/math] and [math]\displaystyle{ \vec{\phi} }[/math]

[math]\displaystyle{ \begin{matrix} \vec{w}=\left[\begin{matrix}{c} w\\lt center\gt \lt math\gt 1pt] \mbox{\Large \lt math\gt \frac{\partial{w}}{\partial{x}} }[/math]} \[math]\displaystyle{ 6pt] \mbox{\Large \lt math\gt \frac{\partial{w}}{\partial{y}} }[/math]} \end{matrix}\right],

\quad \vec{\phi}=\left[\begin{matrix}{c}

\phi \[math]\displaystyle{ 1pt] \mbox{\Large \lt math\gt \frac{\partial{\phi}}{\partial{x}} }[/math]} \[math]\displaystyle{ 6pt] \mbox{\Large \lt math\gt \frac{\partial{\phi}}{\partial{y}} }[/math]} \end{matrix}\right], \end{matrix}</math>

are vectors representing the values of the displacement and potential, respectively, and their horizontal derivatives at the [math]\displaystyle{ q }[/math] nodes and [math]\displaystyle{ \vec{\chi} }[/math] is a vector of the associated basis functions. We can express [math]\displaystyle{ \vec{\chi}^{T}(\mathbf{x}) }[/math] in terms of the FEM basis functions for non-conforming square element denoted by [math]\displaystyle{ \mathbf{N}_{d}\left( \mathbf{x} \right) }[/math], (Wang04 and \cite[pp. 229--243)]{Petyt}) as

[math]\displaystyle{ \vec{\chi}^{T}(\mathbf{x}) = \left( \sum_{d=1}^{N}\mathbf{N}_{d}\left( \mathbf{x} \right) \left[ o\right]_{d}\right) (chi_expansion) }[/math]

where [math]\displaystyle{ N }[/math] is the number of panels and the matrix [math]\displaystyle{ \left[ o\right]_{d} }[/math] is the \emph{assembler matrix \index{Finite Element Method (FEM)!assembler matrix (see also assembler matrix)} \index{assembler matrix}}. Both [math]\displaystyle{ \mathbf{N}_{d}(\mathbf{x}) }[/math] (which is a [math]\displaystyle{ 1 \times 12 }[/math] matrix) and [math]\displaystyle{ \left[ o\right]_{d} }[/math] (which is a [math]\displaystyle{ 12\times 3q }[/math] matrix) are discussed in detail by Wang and Meylan Wang04. Following Wang and Meylan Wang04 we call the square element {\it panel}.

The equation governing the motion of the plate ((maineqn)) and the time-harmonic non-dimensionalised versions of the boundary conditions at the edge of the plate ((boundary1)), ((boundary2)) are equivalent to the variational principle Porter_porter04,AppliedOR01 and \cite[pp. 185--187]{Hildebrand65}, [math]\displaystyle{ \delta L = 0 }[/math] where

[math]\displaystyle{ \begin{matrix} L(w) & = & \frac{1}{2} \int_{\Delta_{0}} \beta \left[ \left(\nabla_h^2 w\right)^2 - 2(1-\nu) \left( \frac{\partial^2 w}{\partial x^2} \frac{\partial^2 w}{\partial y^2} - \left( \frac{\partial^2 w}{\partial x \partial y} \right)^2 \right) \right] \\ & & \qquad + \left( 1- k \gamma \right) w^{2} - 2 k w \phi|_{z=0} \,d\mathbf{x}. \end{matrix} }[/math]

The three terms in the integrand above represent, respectively, the strain energy of the plate, the acceleration, and the dynamic pressure on the plate. Apart from the plate equation itself other natural conditions of [math]\displaystyle{ \delta L = 0 }[/math] are the free edge conditions described by ((boundary1)), and ((boundary2)). Thus, using the variational principle means that the edge conditions are satisfied indirectly as part of the approximation.

If we now substitute the approximation for [math]\displaystyle{ w }[/math] and [math]\displaystyle{ \phi }[/math] in ((expansion)) into the above and minimise (i.e. apply [math]\displaystyle{ \delta L = 0 }[/math]) we obtain

[math]\displaystyle{ \beta \,\mathbb{K}\,\mathbf{\vec{w}}+(1- k \gamma) \,\mathbb{M}\,\mathbf{\vec{w}} = k \,\mathbb{M}\,\mathbf{\vec{\phi}}, (comp_eqn) }[/math]

where the stiffness matrix [math]\displaystyle{ \mathbb{K} }[/math], is given by

[math]\displaystyle{ \begin{matrix} \mathbb{K} & = & \int_{\Delta _{0}}\Bigg[ \frac{\partial ^{2}\vec{\chi}}{\partial x^{2}} \frac{\partial ^{2}\vec{\chi}^{T}}{\partial x^{2}} + \frac{\partial ^{2} \vec{\chi}}{\partial y^{2}} \frac{\partial ^{2}\vec{\chi}^{T}}{\partial y^{2}} \\ & & + 2 \left( 1-\nu \right) \frac{\partial ^{2}\vec{\chi}}{\partial x\partial y} \frac{\partial ^{2}\vec{\chi}^{T}}{\partial x\partial y} +\nu \frac{ \partial ^{2}\vec{\chi}}{\partial x^{2}} \frac{\partial ^{2}\vec{\chi}^{T}}{\partial y^{2}} +\nu \frac{\partial ^{2}\vec{\chi}}{\partial y^{2}} \frac{\partial^{2}\vec{\chi}^{T}}{\partial x^{2}}\Bigg] d\mathbf{x} \end{matrix} }[/math]

and the mass matrix [math]\displaystyle{ \mathbb{M} }[/math] is given by

[math]\displaystyle{ \mathbb{M}=\int_{\Delta _{0}}\vec{\chi}(\mathbf{x})\vec{\chi}^{T}(\mathbf{x} )\,d\mathbf{x}. }[/math]

The integral equation ((bem_eq_2)) is transformed by substituting the approximations for [math]\displaystyle{ w }[/math] and [math]\displaystyle{ \phi }[/math] given by ((expansion)) to obtain

[math]\displaystyle{ \vec{\chi}^{T}(\mathbf{x})\vec{\phi} = \vec{\chi}^{T}(\mathbf{x})\vec{\phi}^{\rm in} + \int_{\Delta _{0}}G_{\mathbf{P}}(\mathbf{x};\mbox{\boldmath\lt math\gt \xi }[/math]})

\left[ k\vec{\chi}^{T}(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]})\vec{\phi}
- \vec{\chi}^{T}(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]})\vec{w}\right]
\,d \mbox{\boldmath[math]\displaystyle{ \xi }[/math]}

</math>

(where where [math]\displaystyle{ \vec{\phi}^{\rm in} }[/math] is the representation of [math]\displaystyle{ {\phi }^{\rm in} }[/math] in the basis functions [math]\displaystyle{ \chi }[/math]) and then multiplying this equation by [math]\displaystyle{ \vec{\chi}(\mathbf{x}) }[/math] and integrating over [math]\displaystyle{ \Delta_{0} }[/math]. This gives us

[math]\displaystyle{ \mathbb{M}\,\mathbf{\vec{\phi}} = \mathbb{M} \vec{\phi}^{\rm in} + k \,\mathbb{G}_{\mathbf{P}}\,\mathbf{\vec{\phi}} - \,\mathbb{G}_{\mathbf{P}}\,\mathbf{\vec{w}}, (phi_plate) }[/math]

where [math]\displaystyle{ \mathbb{G}_{\mathbf{P}} }[/math], is given by

[math]\displaystyle{ \mathbb{G}_{\mathbf{P}} = \int_{\Delta _{0}}\int_{\Delta _{0}}\vec{\chi}(\mathbf{x}) G_{\mathbf{P}}(\mathbf{x},\mbox{\boldmath\lt math\gt \xi }[/math]})

\vec{\chi}^{T}(\mathbf{x})\,d\mathbf{x} d\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}.

</math>

We can use equation ((chi_expansion)) to express the matrices [math]\displaystyle{ \mathbb{K} \lt math\gt , }[/math]\mathbb{M}[math]\displaystyle{ , and }[/math]\mathbb{G}_{\mathbf{P}}</math> as

[math]\displaystyle{ \begin{matrix} \mathbb{K} &=&\sum_{d=1}^{N}\left[ o\right] _{d}^{T}[\int_{\Delta _{d}}\left( \frac{\partial ^{2}\mathbf{N}_{d}^{T}}{\partial x^{2}}\,\frac{ \partial ^{2}\mathbf{N}_{d}}{\partial x^{2}}\right) +\nu \,\left( \frac{ \partial ^{2}\mathbf{N}_{d}^{T}}{\partial x^{2}}\,\frac{\partial ^{2}\mathbf{ N}_{d}}{\partial y^{2}}+\frac{\partial ^{2}\mathbf{N}_{d}^{T}}{\partial y^{2} }\,\frac{\partial ^{2}\mathbf{N}_{d}}{\partial x^{2}}\right) \\ &&+2\left( 1-\nu \right) \,\left( \frac{\partial ^{2}\mathbf{N}_{d}^{T}}{ \partial x\partial y}\,\frac{\partial ^{2}\mathbf{N}_{d}}{\partial x\partial y}\right) +\left( \frac{\partial ^{2}\mathbf{N}_{d}^{T}}{\partial y^{2}}\, \frac{\partial ^{2}\mathbf{N}_{d}}{\partial y^{2}}\right) d\mathbf{x]}\left[ o\right] _{d}, \end{matrix} }[/math] [math]\displaystyle{ \mathbb{M}=\sum_{d=1}^{N}\left[ o\right] _{d}^{T}\left[ \int_{\Delta _{d}} \mathbf{N}_{d}^{T}\left( \mathbf{x}\right) \mathbf{\,N}_{d}\left( \mathbf{x} \right) \,d\mathbf{x}\right] \left[ o\right] _{d}, }[/math]

and

[math]\displaystyle{ \mathbb{G}_{\mathbf{P}}=\sum_{e=1}^{p}\sum_{d=1}^{p}\left[ o\right] _{e}^{T} \left[ \int_{\Delta _{e}}\int_{\Delta _{d}}\mathbf{N}_{e}^{T}\left( \mathbf{x }\right) G_{\mathbf{P}}(\mathbf{x},\mathbf{\xi })\mathbf{N}_{d}\left( \mathbf{\xi }\right) \,d\mathbf{x}d\mathbf{\xi }\right] \left[ o\right] _{d}. }[/math]

These are the equations which are used to find to calculate the matrices [math]\displaystyle{ \mathbb{K} }[/math], [math]\displaystyle{ \mathbb{M} }[/math], and [math]\displaystyle{ \mathbb{G} }[/math]. The solution for [math]\displaystyle{ \vec{w} }[/math] and [math]\displaystyle{ \vec{\phi} }[/math] is then found by solving equation ((comp_eqn)) simultaneously with ((phi_plate)).

Accelerating the Convergence of the Periodic Green Function

The spatial representation of the periodic Green function given by equation ((G_periodic)) is slowly convergent and in the far field the terms decay in magnitude like [math]\displaystyle{ O(n^{-1/2}) }[/math]. In this section we show how to accelerate the convergence. We begin with the asymptotic approximation of the three-dimensional Green function ((green)) far from the source point,

[math]\displaystyle{ G(\mathbf{x},0;\mbox{\boldmath\lt math\gt \xi }[/math]}) \sim -\frac{ik}{2}

\,H_{0}( k |\mathbf{x}-\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}|), 
\qquad \mbox{as [math]\displaystyle{ |\mathbf{x}-\mbox{\boldmath }[/math]\xi[math]\displaystyle{ }| \to \infty }[/math]},
 (asyG)

</math>

Weh_Lait where [math]\displaystyle{ H_0 \equiv H_{0}^{(1)} }[/math] is the Hankel function of the first kind of order zero abr_ste. In Linton Linton98 various methods were described in which the convergence of the periodic Green functions was improved. One such method, which suits the particular problem being considered here, involves writing the periodic Green function as

[math]\displaystyle{ \begin{matrix} G_{\mathbf{P}} (\mathbf{x};\mbox{\boldmath\lt math\gt \xi }[/math]})

&\! = &\!\!\!\! \displaystyle{\sum_{m=-\infty}^{\infty}} \!
\left[ G\left(\mathbf{x};\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}+(0,ml)\right)  
+ \frac{ik}{2} H_{0} \Big(k\sqrt{\left( X+cl\right)^2 + Y_{m}^2}\Big)
\right] \exp^{im\sigma l}
\\
&\! &\! -\displaystyle{\sum_{m=-\infty}^{\infty}}
\frac{ik}{2}H_{0} \Big(k\sqrt{ (X+cl)^2 + Y_{m}^2 }\Big)
\exp^{im\sigma l}
 (near_accelerated)

\end{matrix}</math>

where [math]\displaystyle{ c }[/math] is a numerical smoothing parameter, introduced to avoid the singularity at [math]\displaystyle{ \mathbf{x} = \mbox{\boldmath }[/math]\xi[math]\displaystyle{ } }[/math] in the Hankel function and

[math]\displaystyle{ X = x-\xi,\quad \mbox{and} \quad Y_{m} = (y-\eta)-ml. }[/math]

Furthermore we use the fact that second slowly convergent sum in ((near_accelerated)) can be transformed to

[math]\displaystyle{ (help_accelerated) -\frac{i}{l} \sum_{m=-\infty}^{\infty} \frac{\exp^{ik \mu_{m} |X+c| }\,\exp^{i \sigma_{m} Y_0}}{\mu_{m}} }[/math]

Linton98,Jorgenson90,Singh90 where [math]\displaystyle{ \sigma_{m} = \sigma + 2 m \pi/l }[/math] and

[math]\displaystyle{ \mu_m = \left[ 1-\left(\frac{\sigma_{m}}{k} \right)^{2} \right]^{\frac{1}{2}}, }[/math]

where the positive real or positive imaginary part of the square root is taken. Combining equations ((near_accelerated)) and ((help_accelerated)) we obtain the accelerated version of the periodic Green function

[math]\displaystyle{ \begin{matrix} G_{\mathbf{P}} (\mathbf{x};\mbox{\boldmath\lt math\gt \xi }[/math]})

&\! = &\!\!\!\! \displaystyle{\sum_{m=-\infty}^{\infty}} \!
\left[ G\left(\mathbf{x};\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}+(0,ml)\right)  
+ \frac{ik}{2} H_{0} \Big(k \sqrt{\left( X+cl\right)^2 + Y_{m}^2}\Big)
\right] \exp^{im\sigma l}
\\
&\! &\! -\frac{i}{l}\displaystyle{\sum_{m=-\infty}^{\infty}
\frac{\exp^{ik \mu_{m}|X+cl| }\exp^{i\sigma_{m}Y_0}}{\mu_{m}}}.
 (Gp_fast)

\end{matrix}</math>

The convergence of the two sums in ((Gp_fast)) depends on the value of [math]\displaystyle{ c }[/math]. For small [math]\displaystyle{ c }[/math] the first sum converges rapidly while the second converges slowly. For large [math]\displaystyle{ c }[/math] the second sum converges rapidly while the first converges slowly. The smoothing parameter [math]\displaystyle{ c }[/math] must be carefully chosen to balance these two effects. Of course, the convergence also depends strongly on how close together the points [math]\displaystyle{ \mathbf{x} }[/math] and [math]\displaystyle{ \mbox{\boldmath }[/math]\xi[math]\displaystyle{ } }[/math] are.

Note that some special combinations of wavelength [math]\displaystyle{ \lambda }[/math] and angle of incidence [math]\displaystyle{ \theta }[/math] cause the periodic Green function to diverge ( Scott98,Jorgenson90). This singularity is closely related to the diffracted waves and will be explained shortly.

The scattered waves (modes)

We begin with the accelerated periodic Green function, equation ((Gp_fast)) setting [math]\displaystyle{ c=0 }[/math] and considering the case when [math]\displaystyle{ X }[/math] is large (positive or negative). We also note that for [math]\displaystyle{ m }[/math] sufficiently small or large [math]\displaystyle{ i\mu_m }[/math] will be negative and the corresponding terms will decay. Therefore

[math]\displaystyle{ (G_p_accelerate) G_{\mathbf{P}} (\mathbf{x};\mbox{\boldmath\lt math\gt \xi }[/math]})

\sim - \frac{i}{l} \sum_{m=-M}^{N} \frac{\exp^{ik\mu_{m}|X|}\, \exp^{i\sigma_{m}Y_0}}
{\mu_{m}}, \qquad \mbox{as [math]\displaystyle{ X \to \pm \infty }[/math]}
 (new_G_p_trunc)

</math>

where the integers [math]\displaystyle{ M }[/math] and [math]\displaystyle{ N }[/math] satisfy the following inequalities

[math]\displaystyle{ (M_N) \left. \begin{matrix} [c]{c} \sigma_{-M-1}\lt -k\lt \sigma_{-M},\\ \sigma_{N}\lt k\lt \sigma_{N+1}. \end{matrix} \right\} (cut_off) }[/math]

Equations ((M_N)) can be written as

[math]\displaystyle{ \frac{l}{2\pi}\left(\sigma+k-2\pi \right) \lt M \lt \frac{l}{2\pi}\left( \sigma+k \right), (diffracted_M) }[/math]

and

[math]\displaystyle{ \frac{l}{2\pi}\left( k - \sigma \right) \gt N \gt \frac{l}{2\pi} \left( k-\sigma - 2\pi \right) (diffracted_N) }[/math]

Linton98. It is obvious that [math]\displaystyle{ G_{\mathbf{P}} }[/math] will diverge if [math]\displaystyle{ \sigma_m = \pm k }[/math]; these values correspond to cut-off frequencies which are an expected feature of periodic structures.

The diffracted waves

The diffracted waves are the plane waves which are observed as [math]\displaystyle{ x \to \pm \infty }[/math]. Their amplitude and form are obtained by substituting the limit of the periodic Green function ((new_G_p_trunc)) as [math]\displaystyle{ x\to\pm\infty }[/math] into the boundary integral equation for the potential ((bem_eq_2)). This gives us

[math]\displaystyle{ \lim_{x\to\pm\infty} \phi^{s} ( \mathbf{x},0 ) = - \frac{i}{l} \sum_{m=-M}^{N} \int_{\Delta_0} \frac{\exp^{ik\mu_{m} |X| } \exp^{i\sigma_{m}Y_0}}{\mu_{m}} \left[ k\phi(\mbox{\boldmath\lt math\gt \xi }[/math]},0)

- w(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}) \right] 
d\mbox{\boldmath[math]\displaystyle{ \xi }[/math]},
 (phi_s)

</math>

where [math]\displaystyle{ \phi^{s} = \phi-\phi^{\rm in} }[/math] is the scattered wave which is composed of a finite number of plane waves. For [math]\displaystyle{ x \to -\infty }[/math] the scattered wave is given by

[math]\displaystyle{ \lim_{x\to-\infty}\phi^{s} (\mathbf{x},0) = A_{m}^{-}\,\exp^{ik\mu_{m}x}\exp^{i\sigma_{m}y}, (phi_m_min) }[/math]

where the amplitudes [math]\displaystyle{ A_{m}^{-} }[/math] are identified from ((phi_s)) as

[math]\displaystyle{ A_{m}^{-} = -\frac{i}{\mu_{m}l} \int_{\Delta_0} \exp^{ik\mu_{m}\xi } \exp^{-i\sigma_{m}\eta} \left[ k\phi\left( \mbox{\boldmath\lt math\gt \xi }[/math]}\right)

- w (\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}) \right] 
d\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}.
 (Ad_m)

</math>

Likewise as [math]\displaystyle{ x \to \infty }[/math] the scattered wave is given by

[math]\displaystyle{ \lim_{x\to\infty}\phi^{s} (\mathbf{x},0) = A_{m}^{+} \exp^{-ik\mu_{m}x} \exp^{i\sigma_{m}y}, (phi_m_plus) }[/math]

where [math]\displaystyle{ A_{m}^{+} }[/math] are

[math]\displaystyle{ A_{m}^{+} = -\frac{i}{\mu_{m}l}\int_{\Delta_0} \exp^{-ik\mu_{m}\xi }\exp^{-i\sigma_{m}\eta} \left[ k\phi (\mbox{\boldmath\lt math\gt \xi }[/math]},0)

- w(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}) \right]
d\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}.
 (Ad_p)

</math>

The diffracted waves propagate at various angles with respect to the normal direction of the array. The angles of diffraction, [math]\displaystyle{ \psi_{m}^{\pm} }[/math], are given by

[math]\displaystyle{ \psi_{m}^{\pm} = \tan^{-1}\left( \frac{\sigma_{m}}{\pm k\mu_{m}}\right). (psi_m) }[/math]

Notice that for [math]\displaystyle{ m=0 }[/math] we have

[math]\displaystyle{ \psi_{0}^{\pm}=\pm\theta, (psi_0) }[/math]

where [math]\displaystyle{ \theta }[/math] is the incident angle. This is exactly as expected since we should always have a transmitted wave which travels in the same direction as the incident wave and a reflected wave which travels in the negative incident angle direction.

The fundamental reflected and transmitted waves

We need to be precise when we determine the wave of order zero at [math]\displaystyle{ x\to\infty }[/math] because we have to include the incident wave. There is always at least one set of propagating waves corresponding to [math]\displaystyle{ m=0 }[/math] which correspond to simple reflection and transmission. The coefficient, [math]\displaystyle{ R }[/math], for the fundamental reflected wave for the [math]\displaystyle{ m=0 }[/math] mode is given by

[math]\displaystyle{ R = A_{0}^{-} = -\frac{i}{\mu_{0}l}\int_{\Delta_0}\exp^{ik (\xi\cos\theta -\eta\sin\theta)}\left[ k\phi(\mbox{\boldmath\lt math\gt \xi }[/math]},0)

- w(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]})\right]
d\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}.
 (R)

</math>

The coefficient, [math]\displaystyle{ T }[/math], for the fundamental transmitted wave for the [math]\displaystyle{ m=0 }[/math] mode is given by

[math]\displaystyle{ T = 1 + A_{0}^{+} = 1 - \frac{i}{\mu_{0}l} \int_{\Delta_0} \exp^{-ik(\xi\cos\theta +\eta\sin\theta)}\left[ k\phi(\mbox{\boldmath\lt math\gt \xi }[/math]},0)

- w(\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}) \right]
d\mbox{\boldmath[math]\displaystyle{ \xi }[/math]}.
 (T)

</math>

Conservation of energy

The diffracted wave, taking into account the correction for [math]\displaystyle{ T }[/math], must satisfy the energy flux equation. This simply says that the energy of the incoming wave must be equal to the energy of the outgoing waves. This gives us

[math]\displaystyle{ \cos\theta = \left( |R|^2+|T|^2 \right) \cos\theta + \sum_{\stackrel{\scriptstyle{m=-M} }{m ~ \neq 0}}^{N} \left( |A_{m}^{-}|^2 \cos\psi_{m}^{-} +|A_{m}^{+}|^2 \cos\psi_{m}^{+} \right). (NRGbal) }[/math]

The energy balance equation ((NRGbal)) can be used as an accuracy check on the numerical results.

Results

We tested the convergence of our accelerated version of the Green function and we use [math]\displaystyle{ c = 0.05 }[/math] and [math]\displaystyle{ 44 }[/math] terms in the first sum (the spatial term) and [math]\displaystyle{ 46 }[/math] terms in the second sum (the spectral term) of ((Gp_fast))). These values were determined by a convergence study, the details of which can be found in Wang wangphd04. These values will be used in all our subsequent calculations. We consider four geometries for the plates which are shown in Figure (floes4jfs).

Scattering from a Dock

Aside from the energy balance equation or wide spacing, it is difficult to compare our results to establish their validity. However, there is one case in which we can make comparisons. If we consider the case when we have the dock boundary condition under the plate (so that [math]\displaystyle{ w=0 }[/math]) and the plates are square and joined then the problem reduces to a two dimensional dock problem which is discussed extensively in \cite[Chapter 2]{Linton_mciver01}. To impose the condition of a dock we simply solve equation ((phi_plate)) setting [math]\displaystyle{ \vec{w} }[/math] to zero (we do not require equation ((comp_eqn))), choosing a square plate (geometry 1) and setting the plate separation to [math]\displaystyle{ l=4 }[/math].

Figure (fig_RnT4stiff) shows the reflection and transmission coefficients for a square plate (geometry 1) with the plate separation [math]\displaystyle{ l=4 }[/math] and the dock boundary condition (crosses) and the solution to the two-dimensional dock problem using the method of Linton_mciver01 (solid and dashed lines). As expected the results agree.

Table (table_stiff) shows the values of the coefficient

[math]\displaystyle{ A_{m}^{\pm} }[/math] for a dock of geometry 1 with [math]\displaystyle{ {\lambda=4} }[/math],

[math]\displaystyle{ {l=6} }[/math] and [math]\displaystyle{ {\theta=\pi/6} }[/math]. These results are given to assist in numerical comparisons. \begin{table} \begin{tabular}{ccc} \hline [math]\displaystyle{ m }[/math] & [math]\displaystyle{ A_{m}^{-} }[/math] & [math]\displaystyle{ A_{m}^{+} }[/math] \\ \hline [math]\displaystyle{ -2 }[/math] & [math]\displaystyle{ -0.214-0.042i }[/math] & [math]\displaystyle{ 0.232+0.023i }[/math] \\ [math]\displaystyle{ -1 }[/math] & [math]\displaystyle{ 0.266-0.268i }[/math] & [math]\displaystyle{ -0.185+0.349i }[/math] \\ [math]\displaystyle{ 0 }[/math] & [math]\displaystyle{ 0.631-0.210i }[/math] & [math]\displaystyle{ -0.702-0.141i }[/math] \\ \end{tabular} \caption[]{The coefficients [math]\displaystyle{ A_{m}^{\pm} }[/math] for the case of a dock of geometry 1 with [math]\displaystyle{ {\lambda=4} }[/math], [math]\displaystyle{ {l=6} }[/math] and [math]\displaystyle{ {\theta=\pi/6} }[/math].} (table_stiff) \end{table}

Scattering from Elastic Plates

We begin with a short table of numerical results. Table (flexible_table) is equivalent to Table (table_stiff) except that the plate is now elastic with [math]\displaystyle{ \beta = 0.1 }[/math] and [math]\displaystyle{ \gamma = 0 }[/math]. As expected the reflected energy is less because the waves can propagate under the elastic plates. \begin{table} \begin{tabular}{ccc} \hline [math]\displaystyle{ m }[/math] & [math]\displaystyle{ A_{m}^{-} }[/math] & [math]\displaystyle{ A_{m}^{+} }[/math] \\ \hline [math]\displaystyle{ -2 }[/math] & [math]\displaystyle{ 0.001+0.014i }[/math] & [math]\displaystyle{ -0.040-0.016i }[/math] \\ [math]\displaystyle{ -1 }[/math] & [math]\displaystyle{ -0.016-0.008i }[/math] & [math]\displaystyle{ -0.070-0.099i }[/math] \\ [math]\displaystyle{ 0 }[/math] & [math]\displaystyle{ -0.058-0.072i }[/math] & [math]\displaystyle{ -0.209-0.582i }[/math] \\ \end{tabular} \caption[]{The coefficients [math]\displaystyle{ A_{m}^{\pm} }[/math] for the case of an elastic plate of geometry 1 with [math]\displaystyle{ {\beta=0.1} }[/math], [math]\displaystyle{ {\gamma=0} }[/math], [math]\displaystyle{ {\lambda=4} }[/math], [math]\displaystyle{ {l=6} }[/math] and [math]\displaystyle{ {\theta=\pi/6} }[/math].} (flexible_table) \end{table}

Figures (fig_square_vartheta_beta0p1) and

(fig_triangle_vartheta_beta0p1) show the amplitudes of the diffracted waves due to the array as a function of the incident angle for plates of geometry one and two respectively with [math]\displaystyle{ \beta=0.1 }[/math], [math]\displaystyle{ \gamma=0 }[/math], [math]\displaystyle{ k=\pi/2 }[/math], and [math]\displaystyle{ l=6 }[/math]. There are 3 pairs of diffracted waves (including the reflected-transmitted pair) for any angle. [math]\displaystyle{ G_{\bf P} }[/math] diverges if [math]\displaystyle{ \sigma_n = \pm k }[/math] which for our values of [math]\displaystyle{ l }[/math] and [math]\displaystyle{ k }[/math] means that [math]\displaystyle{ \theta = \pm 0.3398 }[/math]. As [math]\displaystyle{ \theta }[/math] moves across these points one of the diffracted waves disappears (at [math]\displaystyle{ \pm\pi/2 }[/math]) and an other appears (at [math]\displaystyle{ \mp\pi/2 }[/math]). In the plots we have plotted [math]\displaystyle{ A^{\pm}_{-2} }[/math] and [math]\displaystyle{ A^{\pm}_{1} }[/math] with the same line style and also [math]\displaystyle{ A^{\pm}_{-1} }[/math] and [math]\displaystyle{ A^{\pm}_{2} }[/math] since they represent diffracted waves which appear and disappear together. Interestingly the result of doing this is to produce smooth curves for [math]\displaystyle{ -\pi/2\lt \theta\lt \pi/2 }[/math]. Figures (fig_square_vartheta_beta0p1) and (fig_triangle_vartheta_beta0p1) show that there is a very strong dependence on the amount of reflected energy as a function of incident angle with small angles giving the smallest reflection and large angles giving the greatest reflection.

Figures (fig_disp_p_square) to (fig_disp_p_trapezoid) show the real part of the displacement for five plates ([math]\displaystyle{ \Delta_{j} }[/math], [math]\displaystyle{ j=-2,-1,0,1,2 }[/math]) of the array for plates of geometry one to four respectively, with [math]\displaystyle{ \beta=0.1 }[/math], [math]\displaystyle{ \gamma=0 }[/math], and [math]\displaystyle{ l=6 }[/math]. The angle of incidence is [math]\displaystyle{ \theta=\pi/6 }[/math]. We consider two values of the wavenumber, [math]\displaystyle{ k=\pi/2 }[/math] (a) and [math]\displaystyle{ k=\pi/4 }[/math] (b). The complex response of the elastic plates is apparent in these figures as is the coupling between the water and the plate. It is also clear that there is a great deal of difference in the individual behaviour of plates of different geometries. This has practical implications for experiments which might be performed on individual ice floes. For example, from these figures it appears that it will be very difficult to make

measurements of  wave spectra using an accelerometer

deployed on an ice floe if the floe size is comparable to the wavelength.

In Figure (energy_vs_k) we consider the total reflected energy [math]\displaystyle{ E_R }[/math] given by

[math]\displaystyle{ E_R = |R|^2 \cos\theta + \sum_{\stackrel{\scriptstyle{m=-M} }{m ~ \neq 0}}^{N} |A_{m}^{-}|^2 \cos\psi_{m}^{-} }[/math]

as a function of [math]\displaystyle{ k }[/math] for [math]\displaystyle{ l= 6 }[/math], [math]\displaystyle{ \beta=0.1 }[/math], and [math]\displaystyle{ \gamma=0 }[/math] for floes of geometry 1 and 2 and for angles of [math]\displaystyle{ \theta = 0 }[/math], [math]\displaystyle{ \pi/6 }[/math], and [math]\displaystyle{ \pi/3 }[/math]. In Figure (energy_vs_k) we have divided [math]\displaystyle{ E_R }[/math] by [math]\displaystyle{ \cos\theta }[/math] to normalise the curves, since if all the energy is reflected [math]\displaystyle{ E_R/\cos\theta=1 }[/math]. This figure shows the kind of important results which we can produce from our model even with some simple calculations. The results show that for [math]\displaystyle{ k\lt 1 }[/math] there is no scattering of energy at all. While it is to be expected that long wavelength waves will not be strongly scattered this figure gives us a quantitative value for the [math]\displaystyle{ k }[/math] below which there is negligible scattering. As [math]\displaystyle{ k }[/math] increases the reflection increases and this increase is much more marked for waves incident at a large angle. This is to be expected since waves which are normally incident can be expected to pass through the plates even for short wavelengths. However, the effect of angle appears to be much stronger than might be expected on a simple geometric argument (which is also what was found in Figures (fig_square_vartheta_beta0p1) and (fig_triangle_vartheta_beta0p1)). Interestingly, the effect of geometry is much less significant than either [math]\displaystyle{ \theta }[/math] or [math]\displaystyle{ k }[/math]. This implies that for practical purposes it may be sufficient to take one or two representative geometries. This seems surprising given the differences in the responses of plates of different geometries shown in Figures (fig_disp_p_square) to (fig_disp_p_trapezoid). Figure (energy_vs_k) can be regarded as a preliminary figure, showing the kinds of calculations which are required to develop a model for wave scattering in the marginal ice zone from the present results.

Concluding Remarks

Motivated by the problem of modelling wave propagation in the marginal ice zone and by the general problem of scattering by arrays of bodies we have presented a solution to the problem of wave scattering by an infinite array of floating elastic plates. The model is based on the linear theory and assumes that the floe submergence is negligible. For this reason it is applicable to low to moderate wave heights and to floes whose size is large compared to the floe thickness. The solution method is similar to that used to solve for a single plate except that the periodic Green function must be used. The periodic Green function is obtained by summing the free-surface

Green function for all [math]\displaystyle{ y }[/math]. However the periodic Green function is slowly

convergent and therefore a method, which is commonly used in Optics, is devised to accelerate the convergence. The acceleration method involves expressing the infinite sum of the Green function in its far-field form. From the far-field representation of the periodic Green Function, we calculate the diffracted wave far from the array and the cut-off frequencies. We have checked our numerical calculations for energy balance and against the limiting case when the plates are rigid and joined where the solution reduces to that of a rigid dock. We have also presented solutions for a range of elastic plate geometries.

The solution method could be extended to water of finite depth using a similar periodic Green function, but in this case based on the free surface Green function for finite depth water. The same problems of slow convergence would arise and a similar method for convergence of the series would be required. We believe that a method similar to the one presented here could be used to accelerate the convergence of the periodic Green function for water of finite depth. Another extension would be to consider a doubly periodic lattice in which a doubly periodic Green function would be required. This should also be able to be computed quickly by methods similar to the ones presented here.

The most important extension of this work concerns using it to construct a scattering model for the marginal ice zone. Such a model would be based on the solutions presented here but would be a far from trivial extension of it. For a model to be realistic the effects which have been induced by the assumption of periodicity would have to be removed. One method to do this would be to consider random spacings of the ice floes and to average the result over these. This should lead to a scattered wave field over all directions rather than discrete scattering angles as well as a reflected and transmitted wave. These solutions would then have to be combined, again with some kind of averaging and wide spacing approximation, to compute the scattering by multiple rows of ice floes.

\bibliography{/home/groups/seaice/bibdata/mike,/home/groups/seaice/bibdata/others,/home/meylan/Papers/Periodic_Cynthia/ALMOSTALL,/home/meylan/Papers/Periodic_Cynthia/Cynthia}

\pagebreak \begin{figure} [p] \begin{center} \includegraphics[scale = 0.7] {array} \caption{A schematic diagram of the periodic array of floating elastic plates.}

(fig_array)

\end{center} \end{figure}

\begin{figure}[ptb] \begin{center} \includegraphics[scale = 0.7

] {floes4jfs} \caption{Diagram of the four plate geometries for which we will calculate solutions. }

(floes4jfs)

\end{center} \end{figure}

\begin{figure} [p] \begin{center} \includegraphics[scale = 0.7]{stiff_compare} \caption{The reflection coefficient [math]\displaystyle{ R }[/math] (solid line) and the transmission coefficient [math]\displaystyle{ T }[/math] (dashed line) as a function of [math]\displaystyle{ k }[/math] for a two-dimensional dock of length 4 for the incident angles shown. The crosses are the same problem solved using the three-dimensional array code with the dock boundary condition and using plates of geometry 1 with [math]\displaystyle{ l=4 }[/math].}

(fig_RnT4stiff)

\end{center} \end{figure}