Difference between revisions of "Eigenfunction Matching Method for Floating Elastic Plates"

From WikiWaves
Jump to navigationJump to search
 
(36 intermediate revisions by 2 users not shown)
Line 1: Line 1:
=Introduction=
+
{{pages to be deleted}}
 +
 
 +
==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
 
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
+
This work is based on [[Kohout et. al. 2006]]. This is a generalisation of the
assume that the first and last plate are semi-infinite. The presentation here does not
+
[[Eigenfunction Matching Method for a Semi-Infinite Floating Elastic Plate]].
 +
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
 
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
 
open water can be considered by taking the limit as the
Line 12: Line 12:
 
the plate boundary conditions are satisfied as auxiliary equations.
 
the plate boundary conditions are satisfied as auxiliary equations.
  
=Formulation and preliminaries=
+
== Equations ==
  
 
We consider the problem of small-amplitude waves which are incident on a set of floating elastic
 
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
+
plates occupying the entire water surface. The submergence of the plates is considered negligible.  
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>y</math> direction, although we allow the waves to be
 
We assume that the problem is invariant in the <math>y</math> direction, although we allow the waves to be
 
incident from an angle.  
 
incident from an angle.  
 
The set of plates consists of two semi-infinite plates, separated by a region which
 
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.  
 
consists of a finite number of plates with variable properties.  
We note that we can simulate open water by setting the plate
+
We also assume that the plate edges are free to move at
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
 
each boundary, although other boundary conditions could easily be considered using
the methods of solution presented here. A schematic diagram of
+
the methods of solution presented here. We begin with the [[Frequency Domain Problem]] for multiple
the problem is shown in Figure~35.
+
[[Floating Elastic Plate|Floating Elastic Plates]]
 +
in the non-dimensional form of [[Tayler_1986a|Tayler 1986]] ([[Dispersion Relation for a Floating Elastic Plate]])
  
==Assumptions and conditions==
+
<center><math>\begin{matrix}
 
 
 
 
We assume that in the fluid region <math>-\infty < x,y < \infty</math> and
 
<math>-h < 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>\Phi</math>
 
which satisfies Laplace's equation in the fluid region, i.e.
 
<center><math>(3)
 
\nabla^2 \Phi =0,\;\;\;\; \mbox{ for } -h < z \leq 0.
 
</math></center>
 
 
 
We consider only incident waves of a single frequency <math>\omega</math>, and we assume that these waves
 
also have a simple harmonic variation with respect to <math>y</math>. 
 
The velocity potential of the wave can therefore be expressed as [[stoker57,fox_squire94]]:
 
<center><math>\begin{matrix}(4)
 
\Phi(x,y,z,t)= \Re\{\phi(x,z)e^{ik_yy}e^{-i\omega t}\}
 
\end{matrix}</math></center>
 
where <math>\phi</math> is the complex-valued potential, <math>k_y</math> is the wave number in the <math>y</math> direction and <math>\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:
 
<center><math>\begin{matrix}(5)
 
\frac{\partial \Phi}{\partial z} = 0 \;\;\;\; \mbox{ at } z = -h.
 
\end{matrix}</math></center>
 
The corresponding elevation of the plates is defined by
 
<math> \Re\{\eta(x)e^{ik_yy}e^{-i\omega t}\}</math> where, using the linear kinematic condition
 
at the free surface
 
<center><math>\begin{matrix}(6)
 
-i\omega \eta = \frac{\partial \phi}{\partial z} \;\;\;\; \mbox{ at } z = 0,
 
\end{matrix}</math></center>
 
[[billingham_king00]].
 
We assume the <math>\mu</math>th elastic plate has mass density <math>\rho_\mu</math> and thickness <math>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
 
<center><math>\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></center>
 
[[wang_meylan04]] where <math>P</math> is the pressure at the surface, <math>D_\mu</math> is the rigidity constant of the
 
<math>\mu</math>th plate and <math>m_\mu = \rho_\mu d_\mu</math>.
 
The dynamic condition given by the linearised Bernoulli equation applies 
 
<center><math>\begin{matrix}(8)
 
-i\omega \phi + \frac{P}{\rho} + g\eta = 0 \;\;\;\; \mbox { at } z = 0, \;\;\;
 
\end{matrix}</math></center>
 
[[stoker57]], where <math>P</math> is the pressure at the water surface and <math>\rho</math> is the water density.
 
Equating Eq.~\eqref{ElasticPlate} and Eq.~\eqref{LinearizedBournoulli} gives
 
 
 
<center><math>\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></center>
 
 
 
Additional constraints apply at the edges of the elastic plates [[fox_squire94]].
 
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
 
<center><math>\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></center>
 
<center><math>\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></center>
 
where <math>\nu</math> is Poisson's constant and <math>l_\mu</math> and <math>r_\mu</math> represent the left and right edge of the <math>\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>L</math>, and the time variables
 
by a time parameter <math>\sqrt{L/g}</math>.
 
We leave open the choice of length parameter <math>L</math>.
 
The non-dimensional variables, (denoted by an overbar) are
 
<math></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></math>
 
The boundary condition given by Eq.~\eqref{ElasticPlate2} can now be non-dimensionally expressed as
 
<center><math>\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></center>
 
where <math>\beta_\mu = \frac{D_\mu}{\rho_\mu gL^4}</math> is referred to as the stiffness constant and
 
<math>\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>\eta</math> using Eq.~(6), Eqs.~\eqref{eq:u}, \eqref{SeaBed2}, 
 
\eqref{IceEdge1}, \eqref{IceEdge2}, and \eqref{ElasticPlate3}
 
become
 
<center><math>\begin{matrix}(13)
 
 
\left(\frac{\partial^2}{\partial x^2} +  
 
\left(\frac{\partial^2}{\partial x^2} +  
\frac{\partial^2}{\partial z^2} - k_y^2\right) \phi = 0 \;\;\;\; \mbox{ for } -h < z \leq 0,  
+
\frac{\partial^2}{\partial z^2} - k_y^2\right) \phi = 0, \;\;\;\; \mbox{ for } -h < z \leq 0,  
 
\end{matrix}</math></center>
 
\end{matrix}</math></center>
<center><math>\begin{matrix}(14)
+
<center><math>\begin{matrix}
\frac{\partial \phi}{\partial z} = 0 \;\;\;\; \mbox{ at } z = - h,
+
\frac{\partial \phi}{\partial z} = 0, \;\;\;\; \mbox{ at } z = - h,
 
\end{matrix}</math></center>
 
\end{matrix}</math></center>
<center><math>\begin{matrix}(15)
+
<center><math>\begin{matrix}
 
\left( \beta_\mu \left(\frac{\partial^2}{\partial x^2} - k^2_y\right)^2  
 
\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 \;\;\;\;  
+
- \gamma_\mu\alpha + 1\right)\frac{\partial \phi}{\partial z} - \alpha\phi = 0, \;\;\;\;  
 
\mbox{ at } z = 0, \;\;\; l_\mu \leq x \leq r_\mu,
 
\mbox{ at } z = 0, \;\;\; l_\mu \leq x \leq r_\mu,
 
\end{matrix}</math></center>
 
\end{matrix}</math></center>
where <math>\alpha = \omega^2</math> and  
+
where <math>\alpha = \omega^2</math>, <math>\beta_\mu</math> and <math>\gamma_\mu</math>
<center><math>\begin{matrix}(16)
+
and the stiffness and mass constant for the <math>\mu</math>th plate. The conditions
\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,
+
at the edges of the plates are
 +
<center><math>\begin{matrix}
 +
\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></center>
 
\end{matrix}</math></center>
<center><math>\begin{matrix}(17)
+
<center><math>\begin{matrix}
\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.
+
\left(\frac{\partial^2}{\partial x^2} - \nu k^2_y\right)\frac{\partial\phi}{\partial z} = 0,
 +
\;\;\;\;\mbox{ at } z = 0, \;\;\; \mbox{ for } x = l_\mu,r_\mu.
 
\end{matrix}</math></center>
 
\end{matrix}</math></center>
 +
where
 +
<math>l_\mu</math> and <math>r_\mu</math> represent the left and right edge of the <math>\mu</math>th plate as
 +
shown in Figure~35.
  
=Method of solution=
+
==Method of solution==
  
==Eigenfunction expansion==
+
===Eigenfunction expansion===
  
We will solve Eqs.~(13) to (17) using an eigenfunction expansion.
+
We will solve the system of equations using an [[:Category:Eigenfunction Matching Method|Eigenfunction Matching Method]].
This method has been applied in many situations for linear water wave problems,
+
The method was developed by [[Fox and Squire 1994]] for the case of a single plate
and the technique is described in [[linton_mciver01]]. The method was developed by  
+
as the research is described in [[Two-Dimensional Floating Elastic Plate]].
[[fox_squire94]] for the case
 
of the elastic plate boundary condition, and subsequently it has been used by [[barrett_squire96,sahoo01,teng01]].  
 
 
We show here how this method can be extended to the case of an arbitrary number of plates.
 
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
 
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  
 
that extra modes are required in order to solve the higher order boundary  
 
conditions at the plate edges.   
 
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  
 
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  
 
of reflected waves, one of which is propagating but the rest of  
Line 170: Line 75:
 
the transmission and reflection are not expanded at opposite ends of the plate.   
 
the transmission and reflection are not expanded at opposite ends of the plate.   
  
===Separation of variables===
+
====Separation of variables====
 
The potential velocity can be written in terms of an infinite series of separated eigenfunctions under
 
The potential velocity can be written in terms of an infinite series of separated eigenfunctions under
 
each elastic plate, of the form  
 
each elastic plate, of the form  
<math></math>\phi = e^{\kappa_\mu x} \cos(k_\mu(z+h)).<math></math>
+
<math>\phi = e^{\kappa_\mu x} \cos(k_\mu(z+h)).</math>
If we apply the boundary conditions given by Eqs.~(14)
+
If we apply the boundary conditions given  
and (15)
+
we obtain the [[Dispersion Relation for a Floating Elastic Plate]]
we obtain
 
 
<center><math>\begin{matrix}
 
<center><math>\begin{matrix}
k_\mu\tan{(k_\mu h)}= & -\frac{\alpha}{\beta_\mu k_\mu^{4}  + 1 - \alpha\gamma_\mu},
+
k_\mu\tan{(k_\mu h)}= & -\frac{\alpha}{\beta_\mu k_\mu^{4}  + 1 - \alpha\gamma_\mu}  
(18)
 
 
\end{matrix}</math></center>
 
\end{matrix}</math></center>
[[fox_squire94]].
+
Solving for <math>k_\mu</math> gives a pure imaginary root
Solving for <math>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, two complex roots (two complex conjugate paired roots
 
with positive imaginary part in all physical situations), an infinite number of positive real roots  
 
with positive imaginary part in all physical situations), an infinite number of positive real roots  
 
which approach <math>{n\pi}/{h}</math> as <math>n</math> approaches infinity, and also the negative of all  
 
which approach <math>{n\pi}/{h}</math> as <math>n</math> approaches infinity, and also the negative of all  
these roots [[fox_squire94]] . We denote the two complex roots with positive imaginary part  
+
these roots ([[Dispersion Relation for a Floating Elastic Plate]]) . We denote the two complex roots with positive imaginary part  
 
by <math>k_\mu(-2)</math> and <math>k_\mu(-1)</math>, the purely imaginary  
 
by <math>k_\mu(-2)</math> and <math>k_\mu(-1)</math>, the purely imaginary  
 
root with positive imaginary part by <math>k_\mu(0)</math> and the real roots with positive imaginary part
 
root with positive imaginary part by <math>k_\mu(0)</math> and the real roots with positive imaginary part
Line 202: Line 104:
 
boundary and by applying the boundary conditions at the edge of each plate.
 
boundary and by applying the boundary conditions at the edge of each plate.
  
==Expressions for the potential velocity==
+
===Expressions for the potential velocity===
  
 
We now expand the potential under each plate using the separation of variables solution.  
 
We now expand the potential under each plate using the separation of variables solution.  
Line 208: Line 110:
 
at <math>M</math> real roots of the dispersion equation.  
 
at <math>M</math> real roots of the dispersion equation.  
 
The potential <math>\phi</math> can now be expressed as the following sum of eigenfunctions:
 
The potential <math>\phi</math> can now be expressed as the following sum of eigenfunctions:
<center><math>(19)
+
<center><math>
 
\phi \approx \left\{
 
\phi \approx \left\{
\begin{matrix}{ll}
+
\begin{matrix}
 
{
 
{
 
Ie^{\kappa_{1}(0)(x-r_1)}\frac{\cos(k_1(0)(z+h))}{\cos(k_1(0)h)} }+\\
 
Ie^{\kappa_{1}(0)(x-r_1)}\frac{\cos(k_1(0)(z+h))}{\cos(k_1(0)h)} }+\\
Line 236: Line 138:
 
Note that we have divided by  <math>\cos{(kh)}</math>,  
 
Note that we have divided by  <math>\cos{(kh)}</math>,  
 
so that the coefficients are normalised by the
 
so that the coefficients are normalised by the
potential at the free surface rather than at the bottom surface.  
+
potential at the free surface rather than at the bottom surface.
  
 
==Expressions for displacement==
 
==Expressions for displacement==
Line 242: Line 144:
 
<center><math>
 
<center><math>
 
\eta \approx \frac{i}{\omega}\left\{
 
\eta \approx \frac{i}{\omega}\left\{
\begin{matrix}{ll}
+
\begin{matrix}
 
{
 
{
 
Ik_1(0)e^{\kappa_{1}(0)(x-r_1)}\tan{(k_1(0)h)} - } \\
 
Ik_1(0)e^{\kappa_{1}(0)(x-r_1)}\tan{(k_1(0)h)} - } \\
Line 258: Line 160:
 
\end{matrix}\right.
 
\end{matrix}\right.
 
</math></center>
 
</math></center>
 
  
 
==Solving via eigenfunction matching==
 
==Solving via eigenfunction matching==
Line 270: Line 171:
  
 
Taking inner products leads to the following equations
 
Taking inner products leads to the following equations
<center><math>(20)
+
<center><math>
\begin{matrix}{rcl}
+
\begin{matrix}
 
{
 
{
\int_{-h}^0 \phi_\mu(r_\mu,z)\cos \frac{m\pi}{h}(z+h) \, dz } &=&  
+
\int_{-h}^0 \phi_\mu(r_\mu,z)\cos \frac{m\pi}{h}(z+h) \, \mathrm{d}z } &=&  
 
{
 
{
\int_{-h}^0 \phi_{\mu+1}(l_{\mu+1},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) \, \mathrm{d}z }\\  
 
{
 
{
\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}{\partial x}(r_\mu,z) \cos \frac{m\pi}{h}(z+h) \, \mathrm{d}z } &=&  
 
{
 
{
\int_{-h}^0 \frac{\partial\phi_{\mu+1}}{\partial x}(l_{\mu+1},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) \, \mathrm{d}z }
 
\end{matrix}
 
\end{matrix}
 
</math></center>
 
</math></center>
 
where <math>m\in[0,M]</math> and <math>\phi_\mu</math> denotes the potential under the <math>\mu</math>th plate, i.e. the expression
 
where <math>m\in[0,M]</math> and <math>\phi_\mu</math> denotes the potential under the <math>\mu</math>th plate, i.e. the expression
for <math>\phi</math> given by Eq.~(19) valid for <math>l_\mu <x<r_\mu</math>.
+
for <math>\phi</math> valid for <math>l_\mu <x<r_\mu</math>.
 
The remaining equations to be solved are given by the two edge conditions satisfied at both
 
The remaining equations to be solved are given by the two edge conditions satisfied at both
 
edges of each plate
 
edges of each plate
<center><math>(21)
+
<center><math>
\begin{matrix}{rcll}
+
\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} }
 
\left(\frac{\partial^3}{\partial x^3} - (2 - \nu)k_y^2\frac{\partial}{\partial x}\right)\frac{\partial\phi_\mu}{\partial z} }
Line 298: Line 199:
  
 
We will show the explicit form of the linear system of equations which arise  
 
We will show the explicit form of the linear system of equations which arise  
when we solve Eqs.~(20) and (21).  
+
when we solve these equations.  
 
Let <math>{\mathbf T}_\mu</math> be a column vector given by  
 
Let <math>{\mathbf T}_\mu</math> be a column vector given by  
 
<math>\left[T_{\mu}(-2), . . ., T_{\mu}(M)\right]^{{\mathbf T}}</math>  
 
<math>\left[T_{\mu}(-2), . . ., T_{\mu}(M)\right]^{{\mathbf T}}</math>  
 
and <math>{\mathbf R}_\mu</math> be a column vector given by   
 
and <math>{\mathbf R}_\mu</math> be a column vector given by   
<math>\left[R_{\mu}(-2) . . . R_{\mu}(M)\right]^{{\mathbf T}}</math>.\\
+
<math>\left[R_{\mu}(-2) . . . R_{\mu}(M)\right]^{{\mathbf T}}</math>.
  
\noindent The equations which arise from matching at the boundary between the first
+
The equations which arise from matching at the boundary between the first
 
and second plate are
 
and second plate are
<center><math>(22)
+
<center><math>
\begin{matrix}{l}
+
\begin{matrix}
 
I{\mathbf C} + {\mathbf M}^{+}_{R_1} {\mathbf R}_1 ={\mathbf M}^{-}_{T_2}  {\mathbf T}_2
 
I{\mathbf C} + {\mathbf M}^{+}_{R_1} {\mathbf R}_1 ={\mathbf M}^{-}_{T_2}  {\mathbf T}_2
 
+ {\mathbf M}^{-}_{R_2} {\mathbf R}_2,\\
 
+ {\mathbf M}^{-}_{R_2} {\mathbf R}_2,\\
Line 316: Line 217:
 
The equations which arise from matching at the boundary of the <math>\mu</math>th and (<math>\mu+1</math>)th plate
 
The equations which arise from matching at the boundary of the <math>\mu</math>th and (<math>\mu+1</math>)th plate
 
boundary (<math>\mu>1</math>) are  
 
boundary (<math>\mu>1</math>) are  
<center><math>(23)
+
<center><math>
\begin{matrix}{l}
+
\begin{matrix}
 
{\mathbf M}^{+}_{T_\mu} {\mathbf T}_\mu +{\mathbf M}^{+}_{R_\mu} {\mathbf R}_\mu   
 
{\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}^{-}_{T_{\mu+1}} {\mathbf T}_{\mu+1}  
Line 327: Line 228:
 
</math></center>
 
</math></center>
 
The equations which arise from matching at the (<math>\Lambda-1</math>)th and <math>\Lambda</math>th boundary are
 
The equations which arise from matching at the (<math>\Lambda-1</math>)th and <math>\Lambda</math>th boundary are
<center><math>(24)
+
<center><math>
\begin{matrix}{l}
+
\begin{matrix}
 
{\mathbf M}^{+}_{T_{\Lambda-1}} {\mathbf T}_{\Lambda-1}  
 
{\mathbf M}^{+}_{T_{\Lambda-1}} {\mathbf T}_{\Lambda-1}  
 
+ {\mathbf M}^{+}_{R_{\Lambda-1}}  {\mathbf R}_{\Lambda-1}  
 
+ {\mathbf M}^{+}_{R_{\Lambda-1}}  {\mathbf R}_{\Lambda-1}  
Line 337: Line 238:
 
</math></center>
 
</math></center>
  
\noindent where <math>{\mathbf M}^{+}_{T_\mu}</math>, <math>{\mathbf M}^{+}_{R_\mu}</math>,
+
where <math>{\mathbf M}^{+}_{T_\mu}</math>, <math>{\mathbf M}^{+}_{R_\mu}</math>,
<math>{\mathbf M}^{-}_{T_\mu}</math>, and <math>{\mathbf M}^{-}_{R_\mu}</math>are <math>(M+1)</math> by <math>(M+3)</math> matrices given by
+
<math>{\mathbf M}^{-}_{T_\mu}</math>, and <math>{\mathbf M}^{-}_{R_\mu}</math>are <math>(M+1)</math> by <math>(M+3)</math> matrices given by
<center><math>(25)
+
<center><math>
\begin{matrix}{l}
+
\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}^{+}_{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) \, \mathrm{d}z}, \\  
 
{
 
{
{\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}^{+}_{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) \, \mathrm{d}z },\\
 
{
 
{
 
{\mathbf M}^{-}_{T_\mu}(m,n) =  {\mathbf M}^{+}_{R_\mu}(m,n)  }\\
 
{\mathbf M}^{-}_{T_\mu}(m,n) =  {\mathbf M}^{+}_{R_\mu}(m,n)  }\\
Line 352: Line 253:
 
</math></center>
 
</math></center>
 
<math>{\mathbf N}^{+}_{T_\mu}</math>, <math>{\mathbf N}^{+}_{R_\mu}</math>,
 
<math>{\mathbf N}^{+}_{T_\mu}</math>, <math>{\mathbf N}^{+}_{R_\mu}</math>,
<math>{\mathbf N}^{-}_{T_\mu}</math>, and <math>{\mathbf N}^{-}_{R_\mu}</math> are given by
+
<math>{\mathbf N}^{-}_{T_\mu}</math>, and <math>{\mathbf N}^{-}_{R_\mu}</math> are given by
 
<center><math>
 
<center><math>
\begin{matrix}{l}
+
\begin{matrix}
 
{\mathbf N}^{\pm}_{T_\mu}(m,n)= -\kappa_\mu(n){\mathbf M}^{\pm}_{T_\mu}(m,n),\\
 
{\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).
 
{\mathbf N}^{\pm}_{R_\mu}(m,n)= \kappa_\mu(n){\mathbf M}^{\pm}_{R_\mu}(m,n).
Line 360: Line 261:
 
</math></center>
 
</math></center>
 
<math>\mathbf{C}</math> is a <math>(M+1)</math> vector which is given by
 
<math>\mathbf{C}</math> is a <math>(M+1)</math> vector which is given by
<center><math>(26)
+
<center><math>
{\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.
+
{\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)\, \mathrm{d}z.
 
</math></center>
 
</math></center>
  
The integrals in Eqs.~(25) and (26) are each solved analytically. Now, for all but the first and <math>\Lambda</math>th plate, Eq.~(21) becomes
+
The integrals in the above equation are each solved analytically. Now, for all but the first and <math>\Lambda</math>th plate, the edge equation becomes
 
<center><math>
 
<center><math>
\begin{matrix}{l}
+
\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,\\
 
{\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.
Line 376: Line 277:
 
the incident wave. This gives us  
 
the incident wave. This gives us  
 
<center><math>
 
<center><math>
\begin{matrix}{l}
+
\begin{matrix}
 
I \left(
 
I \left(
\begin{matrix}{c}
+
\begin{matrix}
 
{\mathbf E}^{+}_{T_1}(1,0)\\  
 
{\mathbf E}^{+}_{T_1}(1,0)\\  
 
{\mathbf E}^{+}_{T_1}(2,0)
 
{\mathbf E}^{+}_{T_1}(2,0)
Line 388: Line 289:
 
and for the <math>\Lambda</math>th plate we have no reflection so
 
and for the <math>\Lambda</math>th plate we have no reflection so
 
<center><math>
 
<center><math>
\begin{matrix}{l}
+
\begin{matrix}
 
{\mathbf E}^{-}_{T_\mu} {\mathbf T}_\mu = 0.\\
 
{\mathbf E}^{-}_{T_\mu} {\mathbf T}_\mu = 0.\\
 
\end{matrix}
 
\end{matrix}
 
</math></center>
 
</math></center>
  
\noindent <math>{\mathbf E}^{+}_{T_\mu}</math>, <math>{\mathbf E}^{+}_{R_\mu}</math>, <math>{\mathbf E}^{-}_{T_\mu}</math> and <math>{\mathbf E}^{-}_{R_\mu}</math>  are 2 by M+3 matrices given by  
+
<math>{\mathbf E}^{+}_{T_\mu}</math>, <math>{\mathbf E}^{+}_{R_\mu}</math>, <math>{\mathbf E}^{-}_{T_\mu}</math> and <math>{\mathbf E}^{-}_{R_\mu}</math>  are 2 by M+3 matrices given by  
 
<center><math>
 
<center><math>
\begin{matrix}{l}
+
\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)\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}^{+}_{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)}),\\
Line 408: Line 309:
 
</math></center>
 
</math></center>
  
 
+
Now, the matching matrix is a <math>(2M+6)\times(\Lambda-1)</math> by  
\noindent Now, the matching matrix is a <math>(2M+6)\times(\Lambda-1)</math> by  
 
 
<math>(2M+1)\times(\Lambda -1)</math> matrix given by
 
<math>(2M+1)\times(\Lambda -1)</math> matrix given by
 
<center><math>
 
<center><math>
 
{\mathbf M} =  
 
{\mathbf M} =  
\left( \begin{matrix}{cccccccccc}
+
\left( \begin{matrix}
 
{\mathbf M}^{+}_{R_1} & -{\mathbf M}^{-}_{T_2} & -{\mathbf M}^{-}_{R_2} &    0    &    0    &        & 0 & 0 & 0 \\  
 
{\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 \\  
 
{\mathbf N}^{+}_{R_1} & -{\mathbf N}^{-}_{T_2} & -{\mathbf N}^{-}_{R_2} &    0    &    0    &        & 0 & 0 & 0 \\  
Line 424: Line 324:
 
</math></center>
 
</math></center>
  
\noindent the edge matrix is a <math>(2M+6)\times(\Lambda-1)</math> by <math>4(\Lambda-1)</math> matrix given by
+
the edge matrix is a <math>(2M+6)\times(\Lambda-1)</math> by <math>4(\Lambda-1)</math> matrix given by
  
 
<center><math>
 
<center><math>
 
{\mathbf E} =  
 
{\mathbf E} =  
\left( \begin{matrix}{cccccccccc}
+
\left( \begin{matrix}
 
   {\mathbf E}^{+}_{R_1} &    0  &          0        &    0    &    0    &        & 0 & 0 & 0 \\
 
   {\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 \\  
Line 441: Line 341:
 
</math></center>
 
</math></center>
  
\noindent and finally the complete system to be solved is given by
+
and finally the complete system to be solved is given by
 
<center><math>
 
<center><math>
\left( \begin{matrix}{c}
+
\left( \begin{matrix}
 
{\mathbf M}\\
 
{\mathbf M}\\
 
{\mathbf E}\\
 
{\mathbf E}\\
 
\end{matrix} \right)
 
\end{matrix} \right)
 
\times
 
\times
\left( \begin{matrix}{c}
+
\left( \begin{matrix}
 
{\mathbf R}_1\\
 
{\mathbf R}_1\\
 
{\mathbf T}_2\\
 
{\mathbf T}_2\\
Line 461: Line 361:
 
\end{matrix} \right)
 
\end{matrix} \right)
 
=
 
=
\left( \begin{matrix}{c}
+
\left( \begin{matrix}
 
-I{\mathbf C}\\
 
-I{\mathbf C}\\
 
\kappa_{1}(0)I{\mathbf C}\\
 
\kappa_{1}(0)I{\mathbf C}\\
Line 476: Line 376:
 
The final system of equations has size <math>(2M+6)\times (\Lambda - 1)</math> by
 
The final system of equations has size <math>(2M+6)\times (\Lambda - 1)</math> by
 
<math>(2M+6)\times (\Lambda - 1)</math>.  
 
<math>(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.
 
 
==Validation of the solutions==
 
 
 
===Oblique waves through a set of elastic plates with uniform properties===
 
 
The solution method we present can be used to solve many of the simpler problems
 
which have been considered in previous works. For example, we can solve for the problem of a single
 
plate surrounded by water, or for a crack between two semi-infinite plate.
 
We choose here to compare our results with the results of [[porter_evans05]],
 
who solved for the reflection and transmission of flexural-gravity
 
waves propagating obliquely through a set of elastic plates separated by narrow parallel cracks.
 
This is equivalent to our problem if the properties
 
are identical for each elastic plate (the plates are of constant thickness, Young's modulus etc.).
 
We have selected this solution to compare with, because the most challenging aspect of our problem
 
is the fact that we have multiple cracks. We have also compared our solution with a single
 
plate surrounded by water and for the problem of a single crack between two plates.
 
However, we do not present these comparisons here.
 
 
The solution of [[porter_evans05]] expresses the potential <math>\phi</math> in terms of a linear combination of the incident wave and certain source functions located at each of the cracks.
 
Along with satisfying the field and boundary conditions, these source functions satisfy the jump conditions in the displacements and gradients across each crack.
 
 
We will briefly present the solution of [[porter_evans05]] in our notation and non-dimensionalisation.
 
[[porter_evans05]] first define a function <math>\chi(x,z)</math> representing
 
outgoing waves as <math>|x|\rightarrow \infty</math> which satisfies
 
<center><math>
 
\begin{matrix}{rcll}
 
(\nabla^2 - k_y^2)\chi &= &0, &-h<z<0, -\infty<x<\infty,(27)
 
\end{matrix}
 
</math></center>
 
<center><math>
 
\begin{matrix}{rcll}
 
{\frac{\partial\chi}{\partial z} } &=&0, &z=-h, -\infty<x<\infty,(28)
 
\end{matrix}
 
</math></center>
 
<center><math>
 
\begin{matrix}{rcll}
 
{\left( \beta \left(\frac{\partial^2}{\partial x^2} - k^2_y\right)^2 -
 
\gamma\alpha + 1\right)\frac{\partial \chi}{\partial z} - \alpha\chi }
 
&=& \delta(x), &z=0,-\infty<x<\infty,(29)
 
\end{matrix}
 
</math></center>
 
This problem can be solved to give
 
<center><math>
 
\chi(x,z) = -i\sum_{n=-2}^\infty\frac{\sin{(k(n)h)}\cos{(k(n)(z-h))}}{2\alpha C_n}e^{-\kappa(n)|x|},(30)
 
</math></center>
 
where
 
<center><math>
 
C_n=\frac{1}{2}\left(h + \frac{(5\beta k(n)^4 + 1 - \alpha\gamma)\sin^2{(k(n)h)}}{\alpha}\right),
 
</math></center>
 
and <math>k(n)</math> are the solutions of the dispersion Eq.~(18)
 
(remembering that the plate properties are all identical so that there is only
 
a single dispersion equation to solve and we have removed the <math>\mu</math> subscript).
 
 
Consequently, the source functions for a single crack at <math>x</math> = 0 can be defined as
 
<center><math>(31)
 
\begin{matrix}{rcl}
 
\psi_s(x,z)&=& \beta(\chi_{xx}(x,z) - \nu k_y^2\chi(x,z)),\\
 
\psi_a(x,z)&=& \beta(\chi_{xxx}(x,z) - \nu_1 k_y^2\chi_x(x,z)),
 
\end{matrix}
 
</math></center>
 
where <math>\nu_1 = 2-\nu</math>.
 
It can easily be shown that <math>\psi_s</math> is symmetric about <math>x = 0</math> and
 
<math>\psi_a</math> is antisymmetric about <math>x = 0</math>.
 
 
Substituting Eq.~\eqref{eq:chi4} into Eq.~\eqref{eq:psi1} gives
 
<center><math>(32)
 
\begin{matrix}{rcl}
 
\psi_s(x,z)&=&
 
{
 
-\frac{\beta}{\alpha}
 
\sum_{n=-2}^\infty
 
\frac{g_n\cos{(k(n)(z+h))}}{2k_{xn}C_n}e^{\kappa_{n}|x|} },\\
 
\psi_a(x,z)&=&
 
{
 
{\rm sgn}(x) i\frac{\beta}{\alpha}\sum_{n=-2}^\infty
 
\frac{g_n'\cos{(k(n)(z+h))}}{2k_{xn}C_n}e^{\kappa_{n}|x|}},
 
\end{matrix}
 
</math></center>
 
where
 
<center><math>
 
\begin{matrix}{rcl}
 
g_n &=& -i\kappa(n)(-\kappa(n)^2 + \nu k_y^2)(\sin{(k(n)h)},\\
 
g'_n&=& \kappa(n)^2(-\kappa(n)^2 + \nu k_y^2)(\sin{(k(n)h)}.
 
\end{matrix}
 
</math></center>
 
[[porter_evans05]] then express the solution to the problem as a linear combination of the
 
incident wave and pairs of source functions at each crack,
 
<center><math>(33)
 
\phi(x,z) =
 
{ Ie^{-\kappa_{1}(0)(x-r_1)}\frac{\cos(k_1(0)(z+h))}{\cos(k_1(0)h)} }
 
+ \sum_{n=1}^{\Lambda-1}(P_n\psi_s(x-r_n,z) + Q_n\psi_a(x-r_n,z))
 
</math></center>
 
where
 
<math>P_n</math> and <math>Q_n</math> are coefficients to be solved which represent the jump in the gradient
 
and elevation respectively of the plates across the crack <math>x = a_j</math>.
 
The coefficients <math>P_n</math> and <math>Q_n</math> are found by applying the edge conditions Eqs.~(16) and (17) to
 
the <math>z</math> derivative of <math>\phi</math> at <math>z=0</math>.
 
 
The reflection and transmission coefficients, <math>R_1(0)</math> and <math>T_\Lambda(0)</math> can be found from Eq.~\eqref{eq:PorterPhi}
 
by taking the limits as <math>x\rightarrow\pm\infty</math> to obtain:
 
<center><math>
 
\begin{matrix}{rcl}
 
R_1(0) e^{-\kappa_(0)r_1}&=& {- \frac{\beta}{\alpha}\sum_{n=1}^{\Lambda-1}
 
\frac{e^{\kappa(0)r_n}}{2k_0C_0}(g'_0Q_n + ig_0P_j)}\\
 
T_\Lambda(0) e^{\kappa(0)l_\Lambda}&=& 1 + {\frac{\beta}{\alpha}\sum_{n=1}^{\Lambda -1}\frac{e^{-\kappa(0)r_n}}{2k_0C_0}(g'_0Q_n - ig_0P_j)}
 
\end{matrix}
 
</math></center>
 
 
Figure 36 shows a comparison between our results and results calculated using
 
the theory of [[porter_evans05]] for <math>\Lambda=2</math> and <math>\Lambda=4</math> with <math>\beta=0.1</math>, <math>\gamma = 0</math> and <math>h=1</math>. The
 
pluses and circles are the results using our theory, and the solid lines are due to [[porter_evans05]].
 
As can be seen from the figure, the two methods are in close agreement.
 
 
===Solution convergence===
 
We can use the solution of [[porter_evans05]] to investigate
 
the convergence of our solution. Table 1 show <math>|T|</math> for our method
 
and for [[porter_evans05]] as a function of <math>M</math> for <math>\Lambda =2</math> and
 
<math>\Lambda =4</math>, <math>\alpha = 5</math>, <math>\beta = 0.1</math>, <math>\gamma=0</math> and <math>h=1</math>. The rate of
 
convergence of the two solutions is almost identical. The accuracy of
 
two decimal places for <math>M=20</math> is sufficient for most practical calculations.
 
 
===Wave tank experiment===
 
 
The solution method is also validated by comparison to a series of experiments which were performed in
 
a two-dimensional wave tank. These experiments were aimed at simulating wave propagation in the marginal
 
ice zone, and results concerned with determining the dispersion equation are described in [[sakai_hanai02]].
 
The wave tank used for the experiment was 26 m long, 0.8 m wide and 0.6 m deep.
 
The waves were generated using a wave-maker set up at the front of the tank and an active wave
 
absorption system was used at the far end of the tank. Elastic sheets were placed on
 
the surface of the wave tank, with negligible gap. The plates occupied a length of 8 m of the tank
 
and the entire width of the tank. We will compare with the experiments which were performed
 
with one 8 m sheet, two 4 m sheets and four 2 m sheets.
 
The elastic plate was 20 mm thick, the Young's modulus <math>E</math> was approximately  <math>650</math> <math>MPa</math> and
 
the density of the plate was <math>914</math> <math>{\rm kgm}^{-3}</math>.
 
The vertical displacement was measured at 25 different points along the plate using ultrasonic sensors.
 
We assume that Poisson's ratio (<math>\nu</math>) <math>=0.3</math>, <math>g=9.8</math> <math>{\rm ms}^{-2}</math> and the density of water (<math>\rho</math>) <math>=1000</math> <math>{\rm kgm}^{-3}</math>.
 
 
Figure~37 shows the results for a single plate with period <math>T=1.4</math> s for three different amplitudes. The figure
 
plots the absolute value of the displacement as predicted by theory (solid line) with the results measured experimentally
 
(pluses).
 
As well as showing good agreement between measurement and theory, this figure also shows that the experimental
 
amplitudes are within  the linear regime because there is little change in the measured results as the amplitude
 
increases (apart from the uniform linear change). Figures~38
 
to 40 show the results for <math>T =</math> 1 s, 1.2 s and 1.4 s
 
for one, two and three plates respectively. The figures show good agreement, with a trend of increasing
 
agreement as the period increased. The only figure where there is poor agreement is Figure~39
 
for 1.4 s. We are uncertain about the origin of this difference. However, the overall agreement is good, especially considering
 
that we are plotting the amplitude of displacement.  We consider this to be strong confirmation that our model is performing adequately.
 
The cusps apparent in Figure 40 for <math>T=1</math> s, are caused by the plates being so short as to be almost rigid and by having a near zero in displacement. The effect, when plotting the absolute value of displacement, is a cusp.
 
 
==Conclusions==
 
 
We have solved for the linear water wave propagation under a set
 
of floating elastic plates. While the problem was two-dimensional,
 
it does allow the waves to be incident at an angle.
 
The elastic plate properties can be set arbitrarily so that the
 
model can also include regions of open water.  The solution
 
method is based on an eigenfunction matching at the boundaries
 
of the plates. We also impose the free-edge conditions at
 
the plate edges, by deriving fewer equations from matching
 
than there are unknowns. This is done in a very natural way
 
because the eigenfunctions under the plate actually contain
 
extra modes. The method is stable, but computationally
 
demanding for large number of plates.
 
 
We have compared our solution  with
 
one derived by [[porter_evans05]], which applies to the case of uniform
 
plate properties, and we found good agreement.
 
We tested the accuracy of our solution and found it was very
 
similar to that of [[porter_evans05]]. To obtain two-decimal places
 
of accuracy required around twenty modes.
 
We compared our solution
 
method to a series of experiments performed in a two-dimensional wave
 
tank. The agreement with the experiments was fairly good, and we
 
believe we can have a high degree of confidence that our solution is correct
 
within the expected numerical errors.
 
 
\section*{Acknowledgements}
 
 
This research was supported by Marsden grant UOO308
 
from the New Zealand government.
 
 
\bibliographystyle{/home/groups/seaice/Alison/Papers/SetPlates/Formating/elsart-harv}
 
\bibliography{/home/groups/seaice/bibdata/mike,/home/groups/seaice/bibdata/others}
 
 
 
 
\begin{table}[htbp!](34)
 
\caption{<math>|T|</math> for <math>\alpha=5</math> , <math>\beta=0.1</math>, <math>\gamma=0</math>, <math>h=1</math>
 
and the values of <math>\Lambda</math> and <math>M</math> shown,
 
calculated by the present method and by the method in
 
[[porter_evans05]] .}
 
\begin{center}
 
\begin{tabular}{c|c|c|c}
 
<math>\Lambda</math> & <math>M</math> & <math>|T|</math> (present method)& <math>|T|</math> ([[porter_evans05]])\\
 
\hline
 
2  &    5 & 0.72897005265395 & 0.68013661602795 \\
 
    &  10 & 0.73710075717437 & 0.73382189306476 \\
 
    &  20 & 0.73943613533854 & 0.73910099180859 \\
 
    &  50 & 0.74014223492682 & 0.74012279625910 \\
 
    &  100 & 0.74024743508561 & 0.74024507931561 \\
 
    &  150 & 0.74026720286310 & 0.74026651629366 \\
 
\hline
 
4  &    5 & 0.78572228609681 & 0.64049634405062 \\
 
    &  10 & 0.81444198211422 & 0.80423931535963 \\
 
    &  20 & 0.82228249776276 & 0.82126508433661 \\
 
    &  50 & 0.82458694969417 & 0.82452862088603 \\
 
    &  100 & 0.82492540871298 & 0.82491836384358 \\
 
    &  150 & 0.82498871994750 & 0.82498666973497
 
\end{tabular}
 
\end{center}
 
\end{table}
 
 
 
 
\begin{figure}[th]
 
\begin{center}
 
\includegraphics[width=.8\textwidth]{Figures/SetPlatesTheta}
 
\end{center}
 
\caption{A schematic diagram showing the set of floating elastic plates and the
 
coordinate systems used in the solution.
 
The three dimensional region is defined by <math>-\infty < x,y < \infty</math> and
 
<math>-h < z \leq 0</math>.
 
<math>I</math> represents the incident wave.
 
<math>R_\mu</math> and <math>T_\mu</math> represent the reflection and transmission coefficients of the <math>\mu</math>th plate,
 
<math>l_\mu</math> and <math>r_\mu</math> represent the left and right edge of the plate <math>\mu</math> and <math>\Lambda</math> represents the last plate. }
 
(35)
 
\end{figure}
 
 
 
 
\begin{figure}[th]
 
\begin{center}
 
\includegraphics[width=.8\textwidth]{Figures/Compare_SourceFunct}
 
\end{center}
 
\caption{<math>|R_1(0)|</math> (pluses) and <math>|T_\Lambda(0)|</math> (circles)
 
versus <math>\alpha</math>.
 
<math>\beta=0.1</math>, <math>\gamma=0</math> and <math>h=1</math>.
 
Figure (a) represents solutions for two plates with the crack at <math>x=0</math> and with <math>\theta = 0</math>.
 
Figure (b) represents solutions for two plates with the crack at <math>x=0</math> and with <math>\theta = \pi/3</math>.
 
Figure (c) represents solutions for four plates with the cracks at <math>x=0,1,2</math> and with <math>\theta = 0</math>.
 
Figure (d) represents solutions for four plates with the cracks at <math>x=0,1,2</math> and with <math>\theta = \pi/12</math>.}
 
(36)
 
\end{figure}
 
 
 
 
\begin{figure}[th]
 
\begin{center}
 
\includegraphics[width=.8\textwidth]{Figures/PlotModExpPl1AllA}
 
\end{center}
 
\caption{<math>|\eta|</math> from our model and from the experiment (pluses) for a single plate with incident amplitude <math>.84</math> (a), <math>1.61</math> (b) and <math>2.47</math> (c). <math>T = 1.4</math>s.}
 
(37)
 
\end{figure}
 
 
 
 
\begin{figure}[th]
 
\begin{center}
 
\includegraphics[width=.8\textwidth]{Figures/PlotCompModExpPl1}
 
\end{center}
 
\caption{<math>|\eta|</math> from our model and from experiment (pluses) for a single plate for the wave periods 1 s (a), 1.2 s (b) and 1.4 s (c).}
 
(38)
 
\end{figure}
 
 
 
 
\begin{figure}[th]
 
\begin{center}
 
\includegraphics[width=.8\textwidth]{Figures/PlotCompModExpPl2}
 
\end{center}
 
\caption{As Figure~38 except for two plates.}
 
(39)
 
\end{figure}
 
 
 
  
\begin{figure}[th]
 
\begin{center}
 
\includegraphics[width=.8\textwidth]{Figures/PlotCompModExpPl4}
 
\end{center}
 
\caption{As Figure~38 except for four plates.}
 
(40)
 
\end{figure}
 
  
 
[[Category:Floating Elastic Plate]]
 
[[Category:Floating Elastic Plate]]
 +
[[Category:Eigenfunction Matching Method]]

Latest revision as of 00:05, 17 October 2009


Introduction

We show here a solution to the problem of wave propagation under many floating elastic plates of variable properties This work is based on Kohout et. al. 2006. This is a generalisation of the Eigenfunction Matching Method for a Semi-Infinite Floating Elastic Plate. 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.

Equations

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. 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 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. We begin with the Frequency Domain Problem for multiple Floating Elastic Plates in the non-dimensional form of Tayler 1986 (Dispersion Relation for a Floating Elastic Plate)

[math]\displaystyle{ \begin{matrix} \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} \frac{\partial \phi}{\partial z} = 0, \;\;\;\; \mbox{ at } z = - h, \end{matrix} }[/math]
[math]\displaystyle{ \begin{matrix} \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], [math]\displaystyle{ \beta_\mu }[/math] and [math]\displaystyle{ \gamma_\mu }[/math] and the stiffness and mass constant for the [math]\displaystyle{ \mu }[/math]th plate. The conditions at the edges of the plates are

[math]\displaystyle{ \begin{matrix} \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} \left(\frac{\partial^2}{\partial x^2} - \nu k^2_y\right)\frac{\partial\phi}{\partial z} = 0, \;\;\;\;\mbox{ at } z = 0, \;\;\; \mbox{ for } x = l_\mu,r_\mu. \end{matrix} }[/math]

where [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.

Method of solution

Eigenfunction expansion

We will solve the system of equations using an Eigenfunction Matching Method. The method was developed by Fox and Squire 1994 for the case of a single plate as the research is described in Two-Dimensional Floating Elastic Plate. 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 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{ \phi = e^{\kappa_\mu x} \cos(k_\mu(z+h)). }[/math] If we apply the boundary conditions given 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} \end{matrix} }[/math]

Solving for [math]\displaystyle{ k_\mu }[/math] 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 (Dispersion Relation for a Floating Elastic Plate) . 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{ \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} { 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{ \begin{matrix} { \int_{-h}^0 \phi_\mu(r_\mu,z)\cos \frac{m\pi}{h}(z+h) \, \mathrm{d}z } &=& { \int_{-h}^0 \phi_{\mu+1}(l_{\mu+1},z)\cos \frac{m\pi}{h}(z+h) \, \mathrm{d}z }\\ { \int_{-h}^0 \frac{\partial\phi_\mu}{\partial x}(r_\mu,z) \cos \frac{m\pi}{h}(z+h) \, \mathrm{d}z } &=& { \int_{-h}^0 \frac{\partial\phi_{\mu+1}}{\partial x}(l_{\mu+1},z) \cos \frac{m\pi}{h}(z+h) \, \mathrm{d}z } \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] 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{ \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 these equations. 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].

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

[math]\displaystyle{ \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{ \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{ \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]

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{ \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) \, \mathrm{d}z}, \\ { {\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) \, \mathrm{d}z },\\ { {\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{ {\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)\, \mathrm{d}z. }[/math]

The integrals in the above equation are each solved analytically. Now, for all but the first and [math]\displaystyle{ \Lambda }[/math]th plate, the edge equation 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]

[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]

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]

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]

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].