Eigenfunction Matching Method for Floating Elastic Plates

Introduction

The problem of two semi-infinite plates of different properties was considered by Fox and Squire 1994 using eigenfunction expansion. Barrett and Squire 1996 extended the solution of Fox and Squire 1994 to two plates of arbitrary properties. We show here a solution to the problem of wave propagation under many floating elastic plates of variable properties which is a generalisation of these previous works and is based on Kohout et. al. 2006. We assume that the first and last plate are semi-infinite. The presentation here does not allow open water (it could be included but makes the formulation more complicated). In any case open water can be considered by taking the limit as the plate thickness tends to zero. The solution is derived using an extended eigenfunction matching method, in which the plate boundary conditions are satisfied as auxiliary equations.

Formulation and preliminaries

We consider the problem of small-amplitude waves which are incident on a set of floating elastic plates occupying the entire water surface. The submergence of the plates is considered negligible. The extension of the method to submerged plates may be possible by modifying the present formulation but this remains a subject for future research. We assume that the problem is invariant in the [math]\displaystyle{ y }[/math] direction, although we allow the waves to be incident from an angle. The set of plates consists of two semi-infinite plates, separated by a region which consists of a finite number of plates with variable properties. We note that we can simulate open water by setting the plate properties, i.e. thickness, to be small or by introducing an additional formulation. To keep the presentation and the computer code which we have developed as simple as possible, we will not present an additional formulation, and we simply set the plate parameters to be sufficiently small if we require open water for any calculations. We also assume that the plate edges are free to move at each boundary, although other boundary conditions could easily be considered using the methods of solution presented here. A schematic diagram of the problem is shown in Figure~35.

Assumptions and conditions

We assume that in the fluid region [math]\displaystyle{ -\infty \lt x,y \lt \infty }[/math] and [math]\displaystyle{ -h \lt z \leq 0 }[/math], the flow is irrotational and inviscid, so that the fluid velocity can be written as the gradient of a velocity potential [math]\displaystyle{ \Phi }[/math] which satisfies Laplace's equation in the fluid region, i.e.

[math]\displaystyle{ (3) \nabla^2 \Phi =0,\;\;\;\; \mbox{ for } -h \lt z \leq 0. }[/math]

We consider only incident waves of a single frequency [math]\displaystyle{ \omega }[/math], and we assume that these waves also have a simple harmonic variation with respect to [math]\displaystyle{ y }[/math] (consider the Frequency Domain Problem) The velocity potential of the wave can therefore be expressed as

[math]\displaystyle{ \begin{matrix}(4) \Phi(x,y,z,t)= \Re\{\phi(x,z)e^{ik_yy}e^{-i\omega t}\} \end{matrix} }[/math]

where [math]\displaystyle{ \phi }[/math] is the complex-valued potential, [math]\displaystyle{ k_y }[/math] is the wave number in the [math]\displaystyle{ y }[/math] direction and [math]\displaystyle{ \Re }[/math] denotes the real part.

We assume that the seabed is impermeable, and therefore the velocity component normal to the sea floor vanishes. Hence, the velocity potential at the sea floor satisfies:

[math]\displaystyle{ \begin{matrix}(5) \frac{\partial \Phi}{\partial z} = 0 \;\;\;\; \mbox{ at } z = -h. \end{matrix} }[/math]

The corresponding elevation of the plates is defined by [math]\displaystyle{ \Re\{\eta(x)e^{ik_yy}e^{-i\omega t}\} }[/math] where, using the linear kinematic condition at the free surface

[math]\displaystyle{ \begin{matrix}(6) -i\omega \eta = \frac{\partial \phi}{\partial z} \;\;\;\; \mbox{ at } z = 0. \end{matrix} }[/math]

We assume the [math]\displaystyle{ \mu }[/math]th elastic plate has mass density [math]\displaystyle{ \rho_\mu }[/math] and thickness [math]\displaystyle{ d_\mu }[/math]. We assume that the amplitude at the free surface is small relative to the wavelength and that the curvature is small and hence linearity can be applied. The equation of motion for the plate is therefore given by the elastic plate equation

[math]\displaystyle{ \begin{matrix} (7) P = D_\mu\left(\frac{\partial^2}{\partial x^2} - k_y^2\right)^2\eta - \omega^2 m_\mu\eta \;\;\;\; \mbox{ at } z = 0, \;\;\; l_\mu \leq x \leq r_\mu, \end{matrix} }[/math]

where [math]\displaystyle{ P }[/math] is the pressure at the surface, [math]\displaystyle{ D_\mu }[/math] is the rigidity constant of the [math]\displaystyle{ \mu }[/math]th plate and [math]\displaystyle{ m_\mu = \rho_\mu d_\mu }[/math]. The dynamic condition given by the linearised Bernoulli equation applies

[math]\displaystyle{ \begin{matrix}(8) -i\omega \phi + \frac{P}{\rho} + g\eta = 0 \;\;\;\; \mbox { at } z = 0, \;\;\; \end{matrix} }[/math]

where [math]\displaystyle{ P }[/math] is the pressure at the water surface and [math]\displaystyle{ \rho }[/math] is the water density. Equating Eq.~\eqref{ElasticPlate} and Eq.~\eqref{LinearizedBournoulli} gives

[math]\displaystyle{ \begin{matrix}(9) D_\mu \left( \frac{\partial^2}{\partial x^2} - k_y^2 \right)^2\eta -\omega^2 m_\mu\eta -i\omega\rho \phi + \rho g \eta = 0 \;\;\;\; \mbox{ at } z = 0, \;\;\; l_\mu \leq x \leq r_\mu. \end{matrix} }[/math]

Additional constraints apply at the edges of the elastic plates. We assume that the plate edges are free, which implies that the bending moment and the shearing forces at the edges are zero. Therefore the edge boundary conditions can be expressed as

[math]\displaystyle{ \begin{matrix}(10) \left(\frac{\partial^3}{\partial x^3} - (2 - \nu)k^2_y\frac{\partial}{\partial x}\right) \eta= 0 \;\;\;\; \mbox{ at } z = 0 \;\;\; \mbox{ for } x = l_\mu,r_\mu, \end{matrix} }[/math] [math]\displaystyle{ \begin{matrix}(11) \left(\frac{\partial^2}{\partial x^2} - \nu k^2_y\right) \eta = 0\mbox{ for } \;\;\;\; \mbox{ at } z = 0 \;\;\; \mbox{ for } x = l_\mu,r_\mu, \end{matrix} }[/math]

where [math]\displaystyle{ \nu }[/math] is Poisson's constant and [math]\displaystyle{ l_\mu }[/math] and [math]\displaystyle{ r_\mu }[/math] represent the left and right edge of the [math]\displaystyle{ \mu }[/math]th plate as shown in Figure~35.

Non-dimensionalising the variables

It is convenient to reduce the number of constants in the equations by non-dimensionalising. We non-dimensionalise by scaling the spatial variables by a length parameter [math]\displaystyle{ L }[/math], and the time variables by a time parameter [math]\displaystyle{ \sqrt{L/g} }[/math]. We leave open the choice of length parameter [math]\displaystyle{ L }[/math]. The non-dimensional variables, (denoted by an overbar) are [math]\displaystyle{ }[/math] \bar{x} = \frac{x}{L}, \bar{y} = \frac{y}{L}, \bar{z} = \frac{z}{L}, \bar{\eta} = \frac{\eta}{L}, \bar{t} = \frac{t}{\sqrt{g/L}} \mbox{ and } \bar{\phi} = \frac{\phi}{L\sqrt{Lg}}.[math]\displaystyle{ }[/math] The boundary condition given by Eq.~\eqref{ElasticPlate2} can now be non-dimensionally expressed as

[math]\displaystyle{ \begin{matrix}(12) \beta_\mu \left(\frac{\partial^2}{\partial \bar{x}^2} - \bar{k_y}^2 \right)^2\bar{\eta} - \bar{\omega}^2\gamma_\mu\bar{\eta} -i\bar{\omega} \bar{\phi} + \bar{\eta}= 0 \;\;\;\; \mbox{ at } z = 0, \;\;\; \bar{l}_\mu \leq \bar{x} \leq \bar{r}_\mu, \end{matrix} }[/math]

where [math]\displaystyle{ \beta_\mu = \frac{D_\mu}{\rho_\mu gL^4} }[/math] is referred to as the stiffness constant and [math]\displaystyle{ \gamma_\mu = \frac{m_\mu}{\rho L} }[/math] is referred to as the mass constant. From here on in, all equations are expressed non-dimensionally, and for simplicity the overbar will be omitted from the dimensionless variables in what follows.

Final Equations

Eliminating [math]\displaystyle{ \eta }[/math] using Eq.~(6), Eqs.~\eqref{eq:u}, \eqref{SeaBed2}, \eqref{IceEdge1}, \eqref{IceEdge2}, and \eqref{ElasticPlate3} become

[math]\displaystyle{ \begin{matrix}(13) \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial z^2} - k_y^2\right) \phi = 0 \;\;\;\; \mbox{ for } -h \lt z \leq 0, \end{matrix} }[/math] [math]\displaystyle{ \begin{matrix}(14) \frac{\partial \phi}{\partial z} = 0 \;\;\;\; \mbox{ at } z = - h, \end{matrix} }[/math] [math]\displaystyle{ \begin{matrix}(15) \left( \beta_\mu \left(\frac{\partial^2}{\partial x^2} - k^2_y\right)^2 - \gamma_\mu\alpha + 1\right)\frac{\partial \phi}{\partial z} - \alpha\phi = 0 \;\;\;\; \mbox{ at } z = 0, \;\;\; l_\mu \leq x \leq r_\mu, \end{matrix} }[/math]

where [math]\displaystyle{ \alpha = \omega^2 }[/math] and

[math]\displaystyle{ \begin{matrix}(16) \left(\frac{\partial^3}{\partial x^3} - (2 - \nu)k^2_y\frac{\partial}{\partial x}\right) \frac{\partial\phi}{\partial z}= 0 \;\;\;\; \mbox{ at } z = 0 \;\;\; \mbox{ for } x = l_\mu,r_\mu, \end{matrix} }[/math] [math]\displaystyle{ \begin{matrix}(17) \left(\frac{\partial^2}{\partial x^2} - \nu k^2_y\right)\frac{\partial\phi}{\partial z} = 0\mbox{ for } \;\;\;\; \mbox{ at } z = 0 \;\;\; \mbox{ for } x = l_\mu,r_\mu. \end{matrix} }[/math]

Method of solution

Eigenfunction expansion

We will solve Eqs.~(13) to (17) using an Eigenfunction Matching Method. The method was developed by We show here how this method can be extended to the case of an arbitrary number of plates. One of the key features in the eigenfunction expansion method for elastic plates is that extra modes are required in order to solve the higher order boundary conditions at the plate edges. The first and last plates are semi-infinite and the middle plates are finite. The potential velocity of the first plate can be expressed as the summation of an incident wave and of reflected waves, one of which is propagating but the rest of which are evanescent and they decay as [math]\displaystyle{ x }[/math] tends to [math]\displaystyle{ -\infty }[/math]. Similarly the potential under the final plate can be expressed as a sum of transmitting waves, one of which is propagating and the rest of which are evanescent and decay towards [math]\displaystyle{ +\infty }[/math]. The potential under the middle plates can be expressed as the sum of transmitting waves and reflected waves, each of which consists of a propagating wave plus evanescent waves which decay as [math]\displaystyle{ x }[/math] decreases or increases respectively. We could combine these waves in the formulation, but because of the exponential growth (or decay) in the [math]\displaystyle{ x }[/math] direction the solution becomes numerically unstable in some cases if the transmission and reflection are not expanded at opposite ends of the plate.

Separation of variables

The potential velocity can be written in terms of an infinite series of separated eigenfunctions under each elastic plate, of the form [math]\displaystyle{ }[/math]\phi = e^{\kappa_\mu x} \cos(k_\mu(z+h)).[math]\displaystyle{ }[/math] If we apply the boundary conditions given by Eqs.~(14) and (15) we obtain the Dispersion Relation for a Floating Elastic Plate

[math]\displaystyle{ \begin{matrix} k_\mu\tan{(k_\mu h)}= & -\frac{\alpha}{\beta_\mu k_\mu^{4} + 1 - \alpha\gamma_\mu} (18) \end{matrix} }[/math]

Solving for [math]\displaystyle{ k_\mu }[/math], this dispersion Eq.~\eqref{DispersionIce} gives a pure imaginary root with positive imaginary part, two complex roots (two complex conjugate paired roots with positive imaginary part in all physical situations), an infinite number of positive real roots which approach [math]\displaystyle{ {n\pi}/{h} }[/math] as [math]\displaystyle{ n }[/math] approaches infinity, and also the negative of all these roots fox_squire94 . We denote the two complex roots with positive imaginary part by [math]\displaystyle{ k_\mu(-2) }[/math] and [math]\displaystyle{ k_\mu(-1) }[/math], the purely imaginary root with positive imaginary part by [math]\displaystyle{ k_\mu(0) }[/math] and the real roots with positive imaginary part by [math]\displaystyle{ k_\mu(n) }[/math] for [math]\displaystyle{ n }[/math] a positive integer. The imaginary root with positive imaginary part corresponds to a reflected travelling mode propagating along the [math]\displaystyle{ x }[/math] axis. The complex roots with positive imaginary parts correspond to damped reflected travelling modes and the real roots correspond to reflected evanescent modes. In a similar manner, the negative of these correspond to the transmitted travelling, damped and evanescent modes respectively. The coefficient [math]\displaystyle{ \kappa_\mu }[/math] is

[math]\displaystyle{ \kappa_\mu(n) = \sqrt{k_\mu(n)^2 + k_y^2}, }[/math]

where the root with positive real part is chosen or if the real part is negative with negative imaginary part. Note that the solutions of the dispersion equation will be different under plates of different properties, and that the expansion is only valid under a single plate. We will solve for the coefficients in the expansion by matching the potential and its [math]\displaystyle{ x }[/math] derivative at each boundary and by applying the boundary conditions at the edge of each plate.

Expressions for the potential velocity

We now expand the potential under each plate using the separation of variables solution. We always include the two complex and one imaginary root, and truncate the expansion at [math]\displaystyle{ M }[/math] real roots of the dispersion equation. The potential [math]\displaystyle{ \phi }[/math] can now be expressed as the following sum of eigenfunctions:

[math]\displaystyle{ (19) \phi \approx \left\{ \begin{matrix} { Ie^{\kappa_{1}(0)(x-r_1)}\frac{\cos(k_1(0)(z+h))}{\cos(k_1(0)h)} }+\\ { \qquad \qquad \sum_{n=-2}^{M}R_1(n) e^{\kappa_{1}(n)(x-r_1)} \frac{\cos(k_1(n)(z+h))}{\cos(k_1(n)h)} },& \mbox{ for } x \lt r_1,\\ { \sum_{n=-2}^{M}T_{\mu}(n) e^{-\kappa_\mu(n)(x-l_\mu)} \frac{\cos(k_\mu(n)(z+h))}{\cos(k_\mu(n)h)} } + \\ { \qquad \qquad \sum_{n=-2}^{M}R_{\mu}(n) e^{\kappa_\mu(n)(x-r_\mu)} \frac{\cos(k_\mu(n)(z+h))}{\cos(k_\mu(n)h)} }, &\mbox{ for } l_\mu\lt x \lt r_\mu,\\ { \sum_{n=-2}^{M}T_{\Lambda}(n)e^{-\kappa_{\Lambda}(n)(x-l_\Lambda)} \frac{\cos(k_\Lambda(n)(z+h))}{\cos(k_\Lambda(n)h)} }, &\mbox{ for } l_\Lambda\lt x, \end{matrix} \right. }[/math]

where [math]\displaystyle{ I }[/math] is the non-dimensional incident wave amplitude in potential, [math]\displaystyle{ \mu }[/math] is the [math]\displaystyle{ \mu^{th} }[/math] plate, [math]\displaystyle{ \Lambda }[/math] is the last plate, [math]\displaystyle{ r_\mu }[/math] represents the [math]\displaystyle{ x }[/math]-coordinate of the right edge of the [math]\displaystyle{ \mu^{th} }[/math] plate, [math]\displaystyle{ l_\mu }[/math] ([math]\displaystyle{ =r_{\mu-1} }[/math]) represents the [math]\displaystyle{ x }[/math]-coordinate of the left edge of the [math]\displaystyle{ \mu^{th} }[/math] plate, [math]\displaystyle{ R_\mu(n) }[/math] represents the reflected potential coefficient of the [math]\displaystyle{ n^{th} }[/math] mode under the [math]\displaystyle{ \mu^{th} }[/math] plate, and [math]\displaystyle{ T_\mu(n) }[/math] represents the transmitted potential coefficient of the [math]\displaystyle{ n^{th} }[/math] mode under the [math]\displaystyle{ \mu^{th} }[/math] plate. Note that we have divided by [math]\displaystyle{ \cos{(kh)} }[/math], so that the coefficients are normalised by the potential at the free surface rather than at the bottom surface.

Expressions for displacement

The displacement is given by

[math]\displaystyle{ \eta \approx \frac{i}{\omega}\left\{ \begin{matrix}{ll} { Ik_1(0)e^{\kappa_{1}(0)(x-r_1)}\tan{(k_1(0)h)} - } \\ { \qquad \sum_{n=-2}^{M}R_{1}(n)k_1(n)e^{\kappa_{1}(n)(x-r_1)} \tan{(k_1(n)h)} }, & \mbox{ for } x \lt r_1, \\ { -\sum_{n=-2}^{M}T_{\mu}(n)k_\mu(n)e^{-\kappa_\mu(n)(x-l_\mu)}\tan{(k_\mu(n)h)} - }\\ { \qquad \sum_{n=-2}^{M}R_{\mu}(n)k_\mu(n)e^{\kappa_\mu(n)(x-r_\mu)} \tan{(k_\mu(n)h)} },& \mbox{ for } l_\mu\lt x \lt r_\mu,\\ { -\sum_{n=-2}^{M}T_{\Lambda}(n)k_\mu(n)e^{-\kappa_\mu(n)(x-l_\Lambda)}\tan{(k_\mu(n)h)} }, &\mbox{ for } l_\Lambda\lt x. \end{matrix}\right. }[/math]

Solving via eigenfunction matching

To solve for the coefficients, we require as many equations as we have unknowns. We derive the equations from the free edge conditions and from imposing conditions of continuity of the potential and its derivative in the [math]\displaystyle{ x }[/math]-direction at each plate boundary. We impose the latter condition by taking inner products with respect to the orthogonal functions [math]\displaystyle{ \cos \frac{m\pi}{h}(z+h) }[/math], where [math]\displaystyle{ m }[/math] is a natural number. These functions are chosen for the following reasons. The vertical eigenfunctions [math]\displaystyle{ \cos k_\mu(n)(z+h) }[/math] are not orthogonal (they are not even a basis) and could therefore lead to an ill-conditioned system of equations. Furthermore, by choosing [math]\displaystyle{ \cos \frac{m\pi}{h}(z+h) }[/math] we can use the same functions to take the inner products under every plate. Finally, and most importantly, the plate eigenfunctions approach [math]\displaystyle{ \cos{(m\pi/h)(z + h)} }[/math] for large [math]\displaystyle{ m }[/math], so that as we increase the number of modes the matrices become almost diagonal, leading to a very well-conditioned system of equations.

Taking inner products leads to the following equations

[math]\displaystyle{ (20) \begin{matrix} { \int_{-h}^0 \phi_\mu(r_\mu,z)\cos \frac{m\pi}{h}(z+h) \, dz } &=& { \int_{-h}^0 \phi_{\mu+1}(l_{\mu+1},z)\cos \frac{m\pi}{h}(z+h) \, dz }\\ { \int_{-h}^0 \frac{\partial\phi_\mu}{\partial x}(r_\mu,z) \cos \frac{m\pi}{h}(z+h) \, dz } &=& { \int_{-h}^0 \frac{\partial\phi_{\mu+1}}{\partial x}(l_{\mu+1},z) \cos \frac{m\pi}{h}(z+h) \, dz } \end{matrix} }[/math]

where [math]\displaystyle{ m\in[0,M] }[/math] and [math]\displaystyle{ \phi_\mu }[/math] denotes the potential under the [math]\displaystyle{ \mu }[/math]th plate, i.e. the expression for [math]\displaystyle{ \phi }[/math] given by Eq.~(19) valid for [math]\displaystyle{ l_\mu \lt x\lt r_\mu }[/math]. The remaining equations to be solved are given by the two edge conditions satisfied at both edges of each plate

[math]\displaystyle{ (21) \begin{matrix} { \left(\frac{\partial^3}{\partial x^3} - (2 - \nu)k_y^2\frac{\partial}{\partial x}\right)\frac{\partial\phi_\mu}{\partial z} } &=&0, & \mbox{ for } z = 0 \mbox{ and } x = l_\mu,r_\mu,\\ { \left(\frac{\partial^2}{\partial x^2} - \nu k_y^2\right)\frac{\partial\phi_\mu}{\partial z} } &=&0, & \mbox{ for } z = 0 \mbox{ and } x = l_\mu,r_\mu. \end{matrix} }[/math]

We will show the explicit form of the linear system of equations which arise when we solve Eqs.~(20) and (21). Let [math]\displaystyle{ {\mathbf T}_\mu }[/math] be a column vector given by [math]\displaystyle{ \left[T_{\mu}(-2), . . ., T_{\mu}(M)\right]^{{\mathbf T}} }[/math] and [math]\displaystyle{ {\mathbf R}_\mu }[/math] be a column vector given by [math]\displaystyle{ \left[R_{\mu}(-2) . . . R_{\mu}(M)\right]^{{\mathbf T}} }[/math].\\

\noindent The equations which arise from matching at the boundary between the first and second plate are

[math]\displaystyle{ (22) \begin{matrix} I{\mathbf C} + {\mathbf M}^{+}_{R_1} {\mathbf R}_1 ={\mathbf M}^{-}_{T_2} {\mathbf T}_2 + {\mathbf M}^{-}_{R_2} {\mathbf R}_2,\\ -\kappa_1(0)I\mathbf{C} + {\mathbf N}^{+}_{R_1} {\mathbf R}_1 = {\mathbf N}^{-}_{T_2} {\mathbf T}_2 + {\mathbf N}^{-}_{R_2} {\mathbf R}_2. \end{matrix} }[/math]

The equations which arise from matching at the boundary of the [math]\displaystyle{ \mu }[/math]th and ([math]\displaystyle{ \mu+1 }[/math])th plate boundary ([math]\displaystyle{ \mu\gt 1 }[/math]) are

[math]\displaystyle{ (23) \begin{matrix} {\mathbf M}^{+}_{T_\mu} {\mathbf T}_\mu +{\mathbf M}^{+}_{R_\mu} {\mathbf R}_\mu ={\mathbf M}^{-}_{T_{\mu+1}} {\mathbf T}_{\mu+1} + {\mathbf M}^{-}_{ R_{\mu+1}} {\mathbf R}_{\mu+1}, \\ {\mathbf N}^{+}_{T_\mu} {\mathbf T}_\mu + {\mathbf N}^{+}_{R_\mu} {\mathbf R}_\mu ={\mathbf N}^{-}_{T_{\mu+1}} {\mathbf T}_{\mu +1} +{\mathbf N}^{-}_{ R_{\mu +1}} {\mathbf R}_{\mu +1}. \end{matrix} }[/math]

The equations which arise from matching at the ([math]\displaystyle{ \Lambda-1 }[/math])th and [math]\displaystyle{ \Lambda }[/math]th boundary are

[math]\displaystyle{ (24) \begin{matrix} {\mathbf M}^{+}_{T_{\Lambda-1}} {\mathbf T}_{\Lambda-1} + {\mathbf M}^{+}_{R_{\Lambda-1}} {\mathbf R}_{\Lambda-1} = {\mathbf M}^{-}_{T_\Lambda } {\mathbf T}_{\Lambda}, \\ {\mathbf N}^{+}_{T_{\Lambda-1}} {\mathbf T}_{\Lambda-1} + {\mathbf N}^{+}_{R_{\Lambda-1}} {\mathbf R}_{\Lambda-1} = {\mathbf N}^{-}_{T_\Lambda } {\mathbf T}_{\Lambda}, \end{matrix} }[/math]

\noindent where [math]\displaystyle{ {\mathbf M}^{+}_{T_\mu} }[/math], [math]\displaystyle{ {\mathbf M}^{+}_{R_\mu} }[/math],

[math]\displaystyle{ {\mathbf M}^{-}_{T_\mu} }[/math], and [math]\displaystyle{ {\mathbf M}^{-}_{R_\mu} }[/math]are [math]\displaystyle{ (M+1) }[/math] by [math]\displaystyle{ (M+3) }[/math] matrices given by

[math]\displaystyle{ (25) \begin{matrix} { {\mathbf M}^{+}_{T_\mu}(m,n) = \int_{-h}^0 e^{-\kappa_\mu(n) (r_\mu-l_\mu )} \frac{\cos \left(k_{\mu}(n) (z+h)\right)}{\cos \left(k_{\mu}(n) h\right)} \cos \left(\frac{m\pi}{h}(z+h)\right) \, dz}, \\ { {\mathbf M}^{+}_{R_\mu}(m,n) = \int_{-h}^0 \frac{\cos \left(k_{\mu}(n) (z+h)\right)}{\cos \left(k_{\mu}(n) h\right)} \cos \left(\frac{m\pi}{h}(z+h)\right) \, dz },\\ { {\mathbf M}^{-}_{T_\mu}(m,n) = {\mathbf M}^{+}_{R_\mu}(m,n) }\\ { {\mathbf M}^{-}_{R_\mu}(m,n) = {\mathbf M}^{+}_{T_\mu}(m,n). } \end{matrix} }[/math]

[math]\displaystyle{ {\mathbf N}^{+}_{T_\mu} }[/math], [math]\displaystyle{ {\mathbf N}^{+}_{R_\mu} }[/math],

[math]\displaystyle{ {\mathbf N}^{-}_{T_\mu} }[/math], and [math]\displaystyle{ {\mathbf N}^{-}_{R_\mu} }[/math] are given by

[math]\displaystyle{ \begin{matrix} {\mathbf N}^{\pm}_{T_\mu}(m,n)= -\kappa_\mu(n){\mathbf M}^{\pm}_{T_\mu}(m,n),\\ {\mathbf N}^{\pm}_{R_\mu}(m,n)= \kappa_\mu(n){\mathbf M}^{\pm}_{R_\mu}(m,n). \end{matrix} }[/math]

[math]\displaystyle{ \mathbf{C} }[/math] is a [math]\displaystyle{ (M+1) }[/math] vector which is given by

[math]\displaystyle{ (26) {\mathbf C}(m)=\int_{-h}^0 \frac{\cos (k_1(0)(z+h))}{\cos (k_1(0)h)} \cos \left(\frac{m\pi}{h}(z+h)\right)\, dz. }[/math]

The integrals in Eqs.~(25) and (26) are each solved analytically. Now, for all but the first and [math]\displaystyle{ \Lambda }[/math]th plate, Eq.~(21) becomes

[math]\displaystyle{ \begin{matrix} {\mathbf E}^{+}_{T_\mu} {\mathbf T}_\mu + {\mathbf E}^{+}_{R_\mu} {\mathbf R}_\mu = 0,\\ {\mathbf E}^{-}_{T_\mu} {\mathbf T}_\mu + {\mathbf E}^{-}_{R_\mu} {\mathbf R}_\mu = 0. \end{matrix} }[/math]

The first and last plates only require two equations, because each has only one plate edge. The equation for the first plate must be modified to include the effect of the incident wave. This gives us

[math]\displaystyle{ \begin{matrix} I \left( \begin{matrix} {\mathbf E}^{+}_{T_1}(1,0)\\ {\mathbf E}^{+}_{T_1}(2,0) \end{matrix} \right) + {\mathbf E}^{+}_{R_1} {\mathbf R}_1 = 0,\\ \end{matrix} }[/math]

and for the [math]\displaystyle{ \Lambda }[/math]th plate we have no reflection so

[math]\displaystyle{ \begin{matrix} {\mathbf E}^{-}_{T_\mu} {\mathbf T}_\mu = 0.\\ \end{matrix} }[/math]

\noindent [math]\displaystyle{ {\mathbf E}^{+}_{T_\mu} }[/math], [math]\displaystyle{ {\mathbf E}^{+}_{R_\mu} }[/math], [math]\displaystyle{ {\mathbf E}^{-}_{T_\mu} }[/math] and [math]\displaystyle{ {\mathbf E}^{-}_{R_\mu} }[/math] are 2 by M+3 matrices given by

[math]\displaystyle{ \begin{matrix} {\mathbf E}^{-}_{T_\mu}(1,n) = (\kappa_\mu(n)^2 - (2 - \nu)k_y^2)(k_{\mu}(n)\kappa_\mu(n)\tan{(k_{\mu}(n)h)}),\\ {\mathbf E}^{+}_{T_\mu}(1,n) = (\kappa_\mu(n)^2 - (2 - \nu)k_y^2)(k_{\mu}(n)\kappa_\mu(n)e^{-\kappa_\mu(n)(r_\mu - l_\mu)}\tan{(k_{\mu}(n)h)}),\\ {\mathbf E}^{-}_{R_\mu}(1,n) = (\kappa_\mu(n)^2 - (2 - \nu)k_y^2)(-k_{\mu}(n)\kappa_\mu(n)e^{\kappa_\mu(n)(l_\mu - r_\mu)}\tan{(k_{\mu}(n)h)}),\\ {\mathbf E}^{+}_{R_\mu}(1,n) = (\kappa_\mu(n)^2 - (2 - \nu)k_y^2)(-k_{\mu}(n)\kappa_\mu(n)\tan{(k_{\mu}(n)h)}),\\ {\mathbf E}^{-}_{T_\mu}(2,n) = (\kappa_\mu(n)^2 - \nu k_y^2)(-k_{\mu}(n)\tan{(k_{\mu}(n)h)}),\\ {\mathbf E}^{+}_{T_\mu}(2,n) = (\kappa_\mu(n)^2 - \nu k_y^2)(-k_{\mu}(n)e^{-\kappa_\mu(n)(r_\mu - l_\mu)}\tan{(k_{\mu}(n)h)}),\\ {\mathbf E}^{-}_{R_\mu}(2,n) = (\kappa_\mu(n)^2 - \nu k_y^2)(-k_{\mu}(n)e^{\kappa_\mu(n)(l_\mu - r_\mu)}\tan{(k_{\mu}(n)h)}),\\ {\mathbf E}^{+}_{R_\mu}(2,n) = (\kappa_\mu(n)^2 - \nu k_y^2)(-k_{\mu}(n)\tan{(k_{\mu}(n)h)}).\\ \end{matrix} }[/math]

\noindent Now, the matching matrix is a [math]\displaystyle{ (2M+6)\times(\Lambda-1) }[/math] by [math]\displaystyle{ (2M+1)\times(\Lambda -1) }[/math] matrix given by

[math]\displaystyle{ {\mathbf M} = \left( \begin{matrix} {\mathbf M}^{+}_{R_1} & -{\mathbf M}^{-}_{T_2} & -{\mathbf M}^{-}_{R_2} & 0 & 0 & & 0 & 0 & 0 \\ {\mathbf N}^{+}_{R_1} & -{\mathbf N}^{-}_{T_2} & -{\mathbf N}^{-}_{R_2} & 0 & 0 & & 0 & 0 & 0 \\ 0 & {\mathbf M}^{+}_{T_2} & {\mathbf M}^{+}_{R_2} & -{\mathbf M}^{-}_{T_3} & -{\mathbf M}^{-}_{R_3} & & 0 & 0 & 0 \\ 0 & {\mathbf N}^{+}_{T_2} & {\mathbf N}^{+}_{R_2} & -{\mathbf N}^{-}_{T_3} & -{\mathbf N}^{-}_{R_3} & & 0 & 0 & 0 \\ & & \vdots & & & \ddots & \\ 0 & 0 & 0 & 0 & 0 & & {\mathbf M}^{+}_{T_{\Lambda - 1}} & {\mathbf M}^{+}_{R_{\Lambda - 1}} & -{\mathbf M}^{-}_{ T_{\Lambda}} \\ 0 & 0 & 0 & 0 & 0 & & {\mathbf N}^{+}_{T_{\Lambda - 1}} & {\mathbf N}^{+}_{R_{\Lambda - 1}} & -{\mathbf N}^{-}_{T_{\Lambda }} \\ \end{matrix} \right), }[/math]

\noindent the edge matrix is a [math]\displaystyle{ (2M+6)\times(\Lambda-1) }[/math] by [math]\displaystyle{ 4(\Lambda-1) }[/math] matrix given by

[math]\displaystyle{ {\mathbf E} = \left( \begin{matrix} {\mathbf E}^{+}_{R_1} & 0 & 0 & 0 & 0 & & 0 & 0 & 0 \\ 0 & {\mathbf E}^{+}_{T_2} & {\mathbf E}^{+}_{R_2} & 0 & 0 & & 0 & 0 & 0 \\ 0 & {\mathbf E}^{-}_{T_2} & {\mathbf E}^{-}_{R_2} & 0 & 0 & & 0 & 0 & 0 \\ 0 & 0 & 0 & {\mathbf E}^{+}_{T_3} & {\mathbf E}^{+}_{R_3} & & 0 & 0 & 0 \\ 0 & 0 & 0 & {\mathbf E}^{-}_{T_3} & {\mathbf E}^{-}_{R_3} & & 0 & 0 & 0 \\ & & \vdots & & & \ddots & \\ 0 & 0 & 0 & 0 & 0 & & {\mathbf E}^{+}_{T_{\Lambda-1}} & {\mathbf E}^{+}_{R_{\Lambda-1}} & 0 \\ 0 & 0 & 0 & 0 & 0 & & {\mathbf E}^{-}_{T_{\Lambda-1}} & {\mathbf E}^{-}_{R_{\Lambda-1}} & 0 \\ 0 & 0 & 0 & 0 & 0 & & 0 & 0 & {\mathbf E}^{-}_{ T_\Lambda} \end{matrix} \right), }[/math]

\noindent and finally the complete system to be solved is given by

[math]\displaystyle{ \left( \begin{matrix} {\mathbf M}\\ {\mathbf E}\\ \end{matrix} \right) \times \left( \begin{matrix} {\mathbf R}_1\\ {\mathbf T}_2\\ {\mathbf R}_2\\ {\mathbf T}_3\\ {\mathbf R}_3\\ \vdots\\ {\mathbf T}_{\Lambda-1}\\ {\mathbf R}_{\Lambda-1}\\ {\mathbf T}_{\Lambda} \end{matrix} \right) = \left( \begin{matrix} -I{\mathbf C}\\ \kappa_{1}(0)I{\mathbf C}\\ 0\\ \vdots\\ -IE^{+}_{T_1}(1,0)\\ -IE^{+}_{T_1}(2,0)\\ 0\\ \vdots \end{matrix} \right). }[/math]

The final system of equations has size [math]\displaystyle{ (2M+6)\times (\Lambda - 1) }[/math] by [math]\displaystyle{ (2M+6)\times (\Lambda - 1) }[/math]. The method of solution we have derived is relatively simple and leads to large systems of equations when we simulate multiple plates. Our aim is to produce code which is simple to develop and which we have a strong degree of confidence is numerically accurate and error free. We do not want to make any kind of wide-spacing approximations since real ice fields always have some small floes which we want to be able to simulate. We have used our method to solve for up to a hundred plates in simulations of wave propagation in the marginal ice zone.

The system of equations has a large number of zero entries, due to the fact that each plate couples only with its nearest neighbour. It seems likely that a more sophisticated method of solution could be developed, which exploits this structure. We have have been unable to find such a method due to the difficulty of including the free edge conditions.