Meylan 2002b: Difference between revisions

From WikiWaves
Jump to navigationJump to search
No edit summary
No edit summary
Line 10: Line 10:


[[Category:Reference]]
[[Category:Reference]]
= Spectral Solution of Time Dependent Shallow Water Hydroelasticity =
\author{M\ls I\ls C\ls H\ls A\ls E\ls L\ns H.\ns M\ls E\ls Y\ls L\ls A\ls N<math>
^{1}</math>}
?? and in revised form ??
The spectral theory of a thin plate floating on shallow water is derived and
used to solve the time-dependent motion. This theory is based on an energy
inner product in which the evolution operator becomes unitary. Two solution
methods are presented. In the first, the solution is expanded in the
eigenfunctions of a self-adjoint operator, which are the incoming wave
solutions for a single frequency. In the second, the scattering theory of
Lax-Phillips is used. The Lax-Phillips scattering solution is suitable for
calculating only the free motion of the plate. However, it determines the
modes of vibration of the plate-water system. These modes, both oscillate
and decay, are found by a complex search algorithm based contour
integration. As well as an application to modelling floating runways, the
spectral-theory for a floating thin plate on shallow water is a solvable
model for more complicated hydroelastic systems.
==Introduction==
Hydroelasticity is the study of immersed or floating elastic bodies in a
fluid. It has a wide range of applications including very large floating
structures ([[Kashiwagihydro98]]), ships ([[Bishop]]), breakwaters (
[[Stoker]]) and sea ice ([[Squire_Review]]). One of the best studied
hydroelastic models is the linear floating thin plate ([[Newman_deform]],
[[ohmatsuVLFS]], [[Kagemoto97]], and [[Kashiwagihydro98]]) because
it models many physical systems, such as a floating runway or an ice floe (
[[jgrfloecirc]]).
The time-dependence in linear hydroelastic problems is usually removed by
solving for a single frequency which we will refer to as the \emph{single
frequency solution}. The solution is normally found by expanding the elastic
body motion in the free modes of vibration and solving the fluid equations
using a Green function ([[Bishop]]). This is analogous to solving for a
rigid floating body using the six rigid modes ([[Sarp_Isa]]). While
alternative methods have been developed ([[Kashiwagibspline]], \cite
{ohmatsuVLFS}, and [[Kagemoto97]]), these are based on exploiting some
property such as a simple geometry or high relative stiffness.
In contrast to the single frequency solutions, solving time-dependent linear
hydroelastic systems remains a major challenge. [[Kashiwagitime]] and
[[Endotime]] have developed a time-stepping procedure; however, this
method results in error growth in time. Since the problem is linear it is
solvable by a spectral method which eliminates the long-time growth of
errors. Furthermore, such a method provides information about the behaviour
of the solution, such as the decay constant of the motion. However, the
spectral theory for linear hydroelasticity has not been developed. For this
reason, spectral type solutions such as [[Ohmatsutimedep]], based on a
Fourier expansion of the single frequency solution, have only solved the
problem in restricted circumstances.
A floating thin plate on shallow water is considered in this paper. This
problem has been chosen for the following reasons: while the single
frequency solution is straightforward ([[Stoker]]), the time-dependent
solution has never been calculated; recently [[OhkusuNamba]], \cite
{OhkusuISOPE}, and [[Ertekinshallow1999]] used it to model a floating
runway; the spectral-theory developed here is a solvable model for more
complex hydroelastic systems.
The spectral theory for a thin plate on shallow water is based on an inner
product which gives the energy of the plate-water system. With respect to
this inner product the evolution operator becomes unitary. Two different
solution methods are derived from this spectral theory. The first method is
based on an expansion of the solution in eigenfunctions of a self-adjoint
operator. These are the single frequency solutions. The second solution
method is based on the scattering theory of Lax-Phillips ([[laxphilips]]
). It provides the solution in terms of a countable number of modes which
have both an oscillation and a decay. These modes are important as they can
be used to characterise the response of the system. With the exception of
[[Hazard]], they have not been investigated for hydroelastic systems.
==Formulation: A Thin Plate on Shallow Water==
Figure \ref{shallow} shows a schematic diagram of the problem. The plate is
infinite in the <math>y</math> direction, so that only the <math>x</math> and <math>z</math> directions are
considered. The <math>x</math> direction is horizontal, the positive <math>z</math> axis points
vertically up, and the plate covers the region <math>-b\leqslant x\leqslant b.</math>
The water is of uniform depth <math>h</math> which is small<math>\ </math>enough that the water
may be approximated as shallow ([[Stoker]]). The amplitudes are assumed
small enough that the linear theory is appropriate, and the plate is
sufficiently thin that the shallow draft approximation may be made (\cite
{OhkusuISOPE}). The solution could be extended to waves incident at an angle
on a infinite two dimensional plate, as described in [[OhkusuISOPE]], but
to keep the treatment straightforward this is not done here.
The mathematical description of the problem follows from [[Stoker]]. The
kinematic condition is
<center><math>
\partial _{t}\zeta =-h\partial _{x}^{2}\phi ,  (1)
</math></center>
where <math>\phi </math> is the velocity potential of the water (averaged over the
depth) and <math>\zeta </math> is the displacement of the water surface or the plate
(from the shallow draft approximation). The equation derived by equating the
pressure at the free surface is
<center><math>
-\rho g\zeta -\rho \partial _{t}\phi =\left\{
\begin{matrix}{c}
0,\;\;x\notin (-b,b), \\
D\partial _{x}^{4}\zeta +\rho ^{\prime }d\partial _{t}^{2}\zeta ,\;\;x\in
(-b,b),
\end{matrix}
\right.  (2)
</math></center>
where <math>D</math> is the bending rigidity of the plate per unit length, <math>\rho </math> is
the density of water, <math>\rho ^{\prime }</math> is the average density of the plate,
<math>g</math> is the acceleration due to gravity, and <math>d</math> is the thickness of the plate
<math>.</math> At the ends of the plate the free edge boundary conditions
<center><math>
\lim_{x\downarrow -b}\partial _{x}^{2}\zeta =\lim_{x\uparrow b}\partial
_{x}^{2}\zeta =\lim_{x\downarrow -b}\partial _{x}^{3}\zeta =\lim_{x\uparrow
b}\partial _{x}^{3}\zeta =0  (3)
</math></center>
are applied since these are common in offshore engineering applications (
[[OhkusuISOPE]]). However the theory which will be developed applies
equally to any of the energy-conserving edge conditions such as clamped or
pinned and there is no need for the boundary conditions to be symmetric.
Equation (3) gives the following implied boundary
conditions for <math>\phi </math>
<center><math>
\lim_{x\downarrow -b}\partial _{x}^{4}\phi =\lim_{x\uparrow b}\partial
_{x}^{4}\phi =\lim_{x\downarrow -b}\partial _{x}^{5}\phi =\lim_{x\uparrow
b}\partial _{x}^{5}\phi =0  (4)
</math></center>
which will be used subsequently.
Non-dimensional variables are now introduced. The space variables are
non-dimensionalised using the water depth <math>h,</math> and the time variables are
non-dimensionalised using <math>\sqrt{h/g}</math>. The non-dimensional variables are
<center><math>
\bar{x}=\frac{x}{h},\;\bar{t}=t\sqrt{\frac{g}{h}},\;\bar{\zeta}=\frac{\zeta
}{h},\;\mathrm{and}\;\bar{\phi}=\frac{\phi }{h^{2}\sqrt{g/h}}.
</math></center>
In these new variables, (1) and (2) become
<center><math>
\partial _{\bar{t}}\bar{\zeta}=-\partial _{\bar{x}}^{2}\bar{\phi}
(5)
</math></center>
and
<center><math>
-\bar{\zeta}-\partial _{\bar{t}}\bar{\phi}=\left\{
\begin{matrix}{c}
0,\;\;\bar{x}\notin (-\bar{b},\bar{b}), \\
\beta \partial _{\bar{x}}^{4}\bar{\zeta}+\gamma \partial _{\bar{t}}^{2}\bar{
\zeta},\;\;\bar{x}\in (-\bar{b},\bar{b}),
\end{matrix}
\right.  (6)
</math></center>
where <math>\beta </math> and <math>\gamma </math> are
<center><math>
\beta =\frac{D}{\rho gh^{4}}\;\;\mathrm{and\ \ }\gamma =\frac{\rho ^{\prime
}d}{\rho h}.
</math></center>
For clarity the overbar is dropped from now on.
The main change in extending the formulation to water of finite depth is
that the velocity potential will be governed by Laplace's equation. This
makes the solution of the problem much more computationally demanding since
Laplace's equation must be solved by a numerical method, for example the
boundary element method. Furthermore, the extension of the spectral theory,
which will be developed here for shallow water, to water of finite depth is
non-trivial and remains a subject for further work.
===Neglecting the inertia term.===
It can be assumed that <math>\left| \gamma \partial _{t}^{2}\zeta \right| \ll
\left| \zeta \right| </math> for the following reasons ([[OhkusuISOPE]]). If we
consider a mode of the displacement <math>\zeta =ae^{i\omega t}</math> (where <math>a</math> is
the amplitude) then <math>\partial _{t}^{2}\zeta =-\omega ^{2}ae^{i\omega t}.</math>
For each frequency, <math>\omega ,</math> there is a corresponding wavelength <math>\lambda
=2\pi /\omega .</math> In the non-dimensional variables the wave speed and water
depth are both unity. Since the water is shallow the wavelength <math>\lambda \gg
1<math> and thus </math>\omega \ll 1.</math> It follows that any shallow water mode must
satisfy <math>\left| \partial _{t}^{2}\zeta \right| \ll \left| \zeta \right| .</math>
Also, <math>\gamma \ll 1</math> since <math>\rho ^{\prime }<\rho </math> (otherwise the plate
would sink) and <math>d\ll h</math> (otherwise the submergence of the plate would not
be negligible). Therefore,
<center><math>
\left| \gamma \partial _{t}^{2}\zeta \right| \ll \left| \zeta \right| ,
</math></center>
and we assume in what follows that the inertia, <math>\gamma \partial
_{t}^{2}\zeta ,</math> is zero.
It should be noted that the inclusion of the inertia term in the spectral
theory which will be developed is difficult because it introduces a time
dependence in the energy inner product.
==The Energy Inner Product==
While equations (5) and (6) are not
complicated they cannot be solved in a simple manner. It is not possible to
Fourier transform in space because of the spatial discontinuity of the
differential equations. The Weiner-Hopf technique cannot be used because the
discontinuities divide the space into three regions. A Laplace
transformation in time can be applied but this leads to non-trivial
equations involving a spatially discontinuous differential equation subject
to arbitrary initial conditions. However, straightforward solutions can be
derived using spectral theory.
The spectral-theory solution of equations (5) and (\ref
{displacement2}) is based on the spectral theory for a unitary operator
(essentially, an operator is unitary if the adjoint is also the inverse). We
therefore require an inner product in which the evolution operator is
unitary. This inner product, since the system is conservative, is derived
from the energy. The potential and displacement both contribute to this
energy and we combine them in a two component vector, <math>U\left( x,t\right) </math>,
given by
<center><math>
U\left( x,t\right) =\left(
\begin{matrix}{c}
\phi (x,t) \\
i\zeta (x,t)
\end{matrix}
\right) .  (7)
</math></center>
The energy consists of the kinetic energy of the water (<math>\propto \left| \phi
_{t}^{2}\right| <math>), the potential energy of the water (</math>\propto \left| \phi
^{2}\right| </math>), and the energy of the plate. The energy inner product for
the two vectors
<center><math>
U_{1}=\left(
\begin{matrix}{c}
\phi _{1} \\
i\zeta _{1}
\end{matrix}
\right) \;\;\mathrm{and\ \ }U_{2}=\left(
\begin{matrix}{c}
\phi _{2} \\
i\zeta _{2}
\end{matrix}
\right)
</math></center>
is given by
<center><math>
\left\langle U_{1},U_{2}\right\rangle _{\mathcal{H}}=\left\langle \partial
_{x}\phi _{1},\partial _{x}\phi _{2}\right\rangle +\left\langle \left(
1+\beta \left( H\left( x-b\right) -H\left( x+b\right) \right) \partial
_{x}^{4}\right) i\zeta _{1},i\zeta _{2}\right\rangle ,
(8)
</math></center>
where <math>H</math> is the Heaviside function. The subscript <math>\mathcal{H}</math> is used to
denote the special inner product and the angle brackets without the <math>
\mathcal{H}</math> denote the standard inner product, i.e.
<center><math>
\left\langle f\left( x\right) ,g\left( x\right) \right\rangle =\int_{-\infty
}^{\infty }f\left( x\right) g^{\ast }\left( x\right) dx.
</math></center>
We now write (5) and (6) as
<center><math>\begin{matrix}
\frac{1}{i}\partial _{t}U &=&\mathcal{P}U  (9) \\
U\left( x,t\right) _{t=0} &=&U_{0}\left( x\right) =\left(
\begin{matrix}{c}
\phi _{0}(x) \\
i\zeta _{0}(x)
\end{matrix}
\right)  \nonumber
\end{matrix}</math></center>
where the operator <math>\mathcal{P}</math> is
<center><math>
\mathcal{P=}\left(
\begin{matrix}{cc}
0 & 1+\beta \left( H\left( x-b\right) -H\left( x+b\right) \right) \partial
_{x}^{4} \\
-\partial _{x}^{2} & 0
\end{matrix}
\right) .
</math></center>
<math>\mathcal{P}</math> is self-adjoint with respect to the inner product (\ref
{energyinnerprod}) since <math>\mathcal{P}</math> satisfies
<center><math>
\left\langle \mathcal{P}U_{1},U_{2}\right\rangle _{\mathcal{H}}=\left\langle
U_{1},\mathcal{P}U_{2}\right\rangle _{\mathcal{H}}
</math></center>
from integration by parts and the boundary conditions at the end of the
plate (3). We can express the solution to (\ref
{selfadjoint2}) as
<center><math>
U\left( x,t\right) =e^{i\mathcal{P}t}U_{0}  (10)
</math></center>
where <math>e^{i\mathcal{P}t}</math> is a unitary evolution operator.
==The self-adjoint solution method(11)==
In this section, a solution for the time dependent motion of the plate-water
system is developed using the theory of self-adjoint operators. To evaluate
equation (10) we require a method to calculate the evolution
operator <math>e^{i\mathcal{P}t}</math>. This can be accomplished by using the
eigenfunctions of the operator <math>\mathcal{P},</math> which are the single frequency
solutions.
===Finding the eigenfunctions(12)===
Since <math>\mathcal{P}</math> is self-adjoint, the eigenvalues, <math>\lambda ,</math> must be
real and therefore<math>\ </math>the eigenfunctions of <math>\mathcal{P}</math> are oscillatory
exponentials outside the region of water covered by the plate. Furthermore,
since the plate is finite, the spectrum (set of eigenvalues) is the entire
real numbers. As is expected for two-component systems, there are two
eigenfunctions associated with each eigenvalue <math>\lambda </math>. We choose
incoming waves from the left (<math>\Phi ^{>})</math> and the right (<math>\Phi ^{<})</math> of
unit amplitude as a basis for the eigenspace since they are the standard
single frequency solutions. They have the following asymptotics,
<center><math>
\lim_{x\rightarrow -\infty }\Phi ^{>}=\left(
\begin{matrix}{c}
e^{i\lambda x} \\
\lambda e^{i\lambda x}
\end{matrix}
\right) +S_{12}\left(
\begin{matrix}{c}
e^{-i\lambda x} \\
\lambda e^{-i\lambda x}
\end{matrix}
\right) \;\;\mathrm{and\ \ }\lim_{x\rightarrow \infty }\Phi
^{>}=S_{11}\left(
\begin{matrix}{c}
e^{i\lambda x} \\
\lambda e^{i\lambda x}
\end{matrix}
\right)
</math></center>
and
<center><math>
\lim_{t\rightarrow -\infty }\Phi ^{<}=S_{22}\left(
\begin{matrix}{c}
e^{-i\lambda x} \\
\lambda e^{-i\lambda x}
\end{matrix}
\right) \;\;\mathrm{and\ \ }\lim_{x\rightarrow \infty }\Phi ^{>}=\left(
\begin{matrix}{c}
e^{-i\lambda x} \\
\lambda e^{-i\lambda x}
\end{matrix}
\right) +S_{21}\left(
\begin{matrix}{c}
e^{i\lambda x} \\
\lambda e^{i\lambda x}
\end{matrix}
\right) ,
</math></center>
where <math>S_{11}</math>, <math>S_{12,}</math> <math>S_{21},</math> and <math>S_{22}</math> are the reflection and
transmission coefficients (which must be determined). These eigenfunctions,
which are analogous to the Jost solutions of Schrodinger's equation (\cite
{Chadan89}), will be used to calculated the time-dependent solution.
We find the eigenfunction <math>\Phi ^{>}(\lambda ,x)</math> by solving (\ref
{kinematic2}) and (6) in each region. The two components, <math>
\phi ^{>}\left( \lambda ,x\right) <math> and </math>i\zeta ^{>}\left( \lambda ,x\right)
,</math> are given by
<center><math>
\phi ^{>}\left( \lambda ,x\right) =\left\{
\begin{matrix}{c}
e^{-i\lambda x}+S_{11}\left( \lambda \right) e^{i\lambda x},\;\;x<-b, \\
\sum\limits_{j=1}^{6}\alpha _{j}e^{\mu _{j}\left( \lambda \right)
x},\;\;-b<x<b, \\
S_{12}\left( \lambda \right) e^{-i\lambda x},\;\;x>b,
\end{matrix}
\right.  (13)
</math></center>
and
<center><math>
i\zeta ^{>}\left( \lambda ,x\right) =\left\{
\begin{matrix}{c}
\lambda e^{-i\lambda x}+\lambda S_{11}\left( \lambda \right) e^{i\lambda
x},\;\;x<-b, \\
\frac{-1}{\lambda }\sum\limits_{j=1}^{6}\mu _{j}\left( \lambda \right)
^{2}\alpha _{j}e^{\mu _{j}\left( \lambda \right) x},\;\;-b<x<b, \\
\lambda S_{12}\left( \lambda \right) e^{-i\lambda x},\;\;x>b,
\end{matrix}
\right.
</math></center>
where the coefficients <math>\mu _{j}\left( \lambda \right) </math> are the six roots
of the equation
<center><math>
\beta \mu ^{6}+\mu ^{2}+\lambda ^{2}=0.  (14)
</math></center>
The values of <math>S_{11}\left( \lambda \right) ,</math> <math>S_{12}\left( \lambda \right)
,<math> and </math>\alpha _{j}</math> are found from the boundary conditions (4)
and the continuity of <math>\phi </math> and <math>\partial _{x}\phi </math> at <math>x=\pm b.</math>
Therefore, to find the eigenfunction <math>\Phi ^{>}\left( \lambda ,x\right) ,</math>
we solve the 8 by 8 linear system
<center><math>
\mathbf{M}\vec{a}\mathbf{=}\vec{b},  (15)
</math></center>
where <math>\mathbf{M}</math> is the matrix
<center><math>
\mathbf{M}=\left(
\begin{matrix}{cccccccc}
\mu _{1}^{4}e^{-\mu _{1}b} & \mu _{2}^{4}e^{-\mu _{2}b} & \mu
_{3}^{4}e^{-\mu _{3}b} & \mu _{4}^{4}e^{-\mu _{4}b} & \mu _{5}^{4}e^{-\mu
_{5}b} & \mu _{6}^{4}e^{-\mu _{6}b} & 0 & 0 \\
\mu _{1}^{5}e^{-\mu _{1}b} & \mu _{2}^{5}e^{-\mu _{2}b} & \mu
_{3}^{5}e^{-\mu _{3}b} & \mu _{4}^{5}e^{-\mu _{4}b} & \mu _{5}^{5}e^{-\mu
_{5}b} & \mu _{6}^{5}e^{-\mu _{6}b} & 0 & 0 \\
\mu _{1}^{4}e^{\mu _{1}b} & \mu _{2}^{4}e^{\mu _{2}b} & \mu _{3}^{4}e^{\mu
_{3}b} & \mu _{4}^{4}e^{\mu _{4}b} & \mu _{5}^{4}e^{\mu _{5}b} & \mu
_{6}^{4}e^{\mu _{6}b} & 0 & 0 \\
\mu _{1}^{5}e^{\mu _{1}b} & \mu _{2}^{5}e^{\mu _{2}b} & \mu _{3}^{5}e^{\mu
_{3}b} & \mu _{4}^{5}e^{\mu _{4}b} & \mu _{5}^{5}e^{\mu _{5}b} & \mu
_{6}^{5}e^{\mu _{6}b} & 0 & 0 \\
e^{-\mu _{1}b} & e^{-\mu _{2}b} & e^{-\mu _{3}b} & e^{-\mu _{4}b} & e^{-\mu
_{5}b} & e^{-\mu _{6}b} & -e^{-i\lambda b} & 0 \\
\mu _{1}e^{-\mu _{1}b} & \mu _{2}e^{-\mu _{2}b} & \mu _{3}e^{-\mu _{3}b} &
\mu _{4}e^{-\mu _{4}b} & \mu _{5}e^{-\mu _{5}b} & \mu _{6}e^{-\mu _{6}b} &
-i\lambda e^{-i\lambda b} & 0 \\
e^{\mu _{1}b} & e^{\mu _{2}b} & e^{\mu _{3}b} & e^{\mu _{4}b} & e^{\mu _{5}b}
& e^{\mu _{6}b} & 0 & -e^{-i\lambda b} \\
\mu _{1}e^{\mu _{1}b} & \mu _{2}e^{\mu _{2}b} & \mu _{3}e^{\mu _{3}b} & \mu
_{4}e^{\mu _{4}b} & \mu _{5}e^{\mu _{5}b} & \mu _{6}e^{\mu _{6}b} & 0 &
i\lambda e^{-i\lambda b}
\end{matrix}
\right)
</math></center>
and <math>\vec{a}</math> and <math>\vec{b}</math> are given by
<center><math>
\vec{a}=\left(
\begin{matrix}{c}
\alpha _{1} \\
\alpha _{2} \\
\alpha _{3} \\
\alpha _{4} \\
\alpha _{5} \\
\alpha _{6} \\
S_{11} \\
S_{12}
\end{matrix}
\right) ,\;\;\;\;\;\vec{b}=\left(
\begin{matrix}{c}
0 \\
0 \\
0 \\
0 \\
e^{i\lambda b} \\
-i\lambda e^{i\lambda b} \\
0 \\
0
\end{matrix}
\right) .
</math></center>
Note that the coefficients <math>S_{11}</math> and <math>S_{12}</math> are contained in <math>\vec{a}.</math>
The eigenfunctions for the wave propagating from the right <math>\Phi ^{<}\left(
\lambda ,x\right) <math> are found similarly. Since </math>S_{11}</math> represents the
amplitude of the reflected wave and <math>S_{12}</math> represents the amplitude of the
transmitted wave, conservation of energy requires that <math>\left| S_{11}\right|
^{2}+\left| S_{12}\right| ^{2}=1.</math> Similarly, since the boundary conditions
are symmetric <math>S_{22}\left( \lambda \right) =S_{11}\left( \lambda \right) </math>
and <math>S_{12}\left( \lambda \right) =S_{21}\left( \lambda \right) .</math>
===Solution with the eigenfunctions===
Equation (10) can be solved by a generalised Fourier transform
based on the eigenfunctions of the operator <math>\mathcal{P}</math>. The
eigenfunctions are orthogonal since <math>\mathcal{P}</math> is self-adjoint, but they
must be normalised. This is accomplished by using the following identity
<center><math>
\int_{0}^{\infty }e^{i\left( \lambda _{1}-\lambda _{2}\right) t}dt=\pi
\delta \left( \lambda _{1}-\lambda _{2}\right)
</math></center>
where <math>\delta </math> is the Dirac delta function. Therefore
<center><math>\begin{matrix}
\left\langle \Phi ^{>}\left( x,\lambda _{1}\right) ,\Phi ^{>}\left(
x,\lambda _{2}\right) \right\rangle _{\mathcal{H}} &=&\pi \delta \left(
\lambda _{1}-\lambda _{2}\right) \lambda _{1}^{2}\left( 1+\left|
S_{11}\right| ^{2}+\left| S_{12}\right| ^{2}\right) \\
&&+\pi \delta \left( \lambda _{1}-\lambda _{2}\right) \lambda _{1}\lambda
_{2}\left( 1+\left| S_{11}\right| ^{2}+\left| S_{12}\right| ^{2}\right) \\
&=&4\pi \delta \left( \lambda _{1}-\lambda _{2}\right) \lambda _{1}^{2},
\nonumber
\end{matrix}</math></center>
using the condition <math>\left| S_{11}\right| ^{2}+\left| S_{12}\right| ^{2}=1.</math>
Similarly
<center><math>
\left\langle \Phi ^{<}\left( x,\lambda _{1}\right) ,\Phi ^{<}\left(
x,\lambda _{2}\right) \right\rangle _{\mathcal{H}}=4\pi \delta \left(
\lambda _{1}-\lambda _{2}\right) \lambda _{1}^{2}
</math></center>
and
<center><math>\begin{matrix}
\left\langle \Phi ^{>}\left( x,\lambda _{1}\right) ,\Phi ^{<}\left(
x,\lambda _{2}\right) \right\rangle _{\mathcal{H}} &=&2\pi \delta \left(
\lambda _{1}-\lambda _{2}\right) \lambda _{1}^{2}\left( S_{11}S_{21}^{\ast
}+S_{12}S_{22}^{\ast }\right) \\
&=&0.  \nonumber
\end{matrix}</math></center>
The generalised Fourier transform which solves the evolution equation (\ref
{unitary}) is
<center><math>\begin{matrix}
U\left( x,t\right) &=&\int_{-\infty }^{\infty }\left\langle U_{0}\left(
x\right) ,\frac{\Phi ^{>}\left( x,\lambda \right) }{4\pi \lambda ^{2}}
\right\rangle _{\mathcal{H}}\Phi ^{>}\left( x,\lambda \right) e^{i\lambda
t}d\lambda  (16) \\
&&+\int_{-\infty }^{\infty }\left\langle U_{0}\left( x\right) ,\frac{\Phi
^{<}\left( x,\lambda \right) }{4\pi \lambda ^{2}}\right\rangle _{\mathcal{H}
}\Phi ^{<}\left( x,\lambda \right) e^{i\lambda t}d\lambda .  \nonumber
\end{matrix}</math></center>
Equation (16) is the cornerstone of the approach. The
integral in equation (16) can be calculated by the fast
Fourier transform while the inner product can calculated by the fast Fourier
transform if the initial condition <math>U_{0}</math> is zero underneath the plate (<math>
-b<x<b).</math>
===Numerical Calculations===
The intention of this paper is to develop the solution methods rather than
describe the physics of the motion and therefore only a few solutions are
presented. From [[OhkusuISOPE]] we assume the plate stiffness is <math>\beta
=2\times 10^{4}<math> and the plate length is </math>b=50</math> throughout. These values are
typical for a floating runway. Figures \ref{incomingspecplot1} and \ref
{incomingspecplot2} show the displacement and potential, respectively, for a
pulse travelling to the right at the times <math>t=0,</math> 30, 60, 90, 120, 150, 180,
210, and 240. The incoming wave pulse was chosen to be a Gaussian in
potential centered at <math>x=-125</math> and sufficiently sharp to be negligible under
the plate,
<center><math>
U_{0}\left( x\right) =\left(
\begin{matrix}{c}
\phi \left( x\right) \\
i\phi ^{\prime }\left( x\right)
\end{matrix}
\right)
</math></center>
where
<center><math>
\phi \left( x\right) =\left\{
\begin{matrix}{c}
e^{-\tfrac{(x+125)^{2}}{350}},\;\;x<-50, \\
0,\;\;x>-50.
\end{matrix}
\right.
</math></center>
At <math>t=0</math> the plate is initially at rest and the wave is to the left of the
plate propagating towards it. From <math>t=30</math> the wave has reached the plate and
the plate begins to undergo a complex bending motion in response to the
incoming wave. The response of the plate in turn induces waves in the
surrounding water which propagate away from the plate to the left and right.
The final picture, <math>t=240</math>, shows the plate at rest with waves now
propagating away from it. The majority of the wave energy has passed under
the plate and continues to propagate to the right. However, the shape of the
outgoing wave profile is markedly different from the incoming wave profile.
Also, there is a significant reflected wave propagating away from the plate
to the left.
Figures \ref{spectral1} and \ref{spectral2} show the evolution of the plate
from an initial displacement in the absence of wave forcing for the times <math>
t=0,</math> 20, 40, 60, 80, 100, 120, 140, and 160. Only the plate displacement is
initially non-zero so that
<center><math>
U_{0}\left( x\right) =\left(
\begin{matrix}{c}
0 \\
i\zeta \left( x\right)
\end{matrix}
\right) .
</math></center>
Figure \ref{spectral1} shows the motion for the symmetric initial plate
displacement
<center><math>
\zeta \left( x\right) =e^{-\tfrac{x^{2}}{350}}.
</math></center>
As the plate evolves the plate vibrates, straightens, and the amplitude
decays. The decay is due to the radiation of energy by the waves generated
in the surrounding water. A complex wave train is produced by the plate
motion and can be seen propagating away from the plate. Figure \ref
{spectral2} shows the motion for the non-symmetric initial plate
displacement
<center><math>
\zeta \left( x\right) =e^{-\tfrac{\left( x-50\right) ^{2}}{350}}.
</math></center>
Again as the plate evolves it straightens, vibrates, and decays and induces
waves in the surrounding water.
==The Lax-Phillips Scattering Solution Method.==
In this section, a solution to the time-dependent motion of the plate-water
system is developed using the Lax-Phillips scattering theory (\cite
{laxphilips}). This solution method will only solve for an initial condition
which is zero outside the plate, i.e. <math>U_{0}\left( x\right) =0</math> if <math>\left|
x\right| >b</math>. However, it calculates the solution by an expansion in a
countable number of modes.
===Lax-Phillips Scattering(17)===
The Lax-Phillips scattering theory will be briefly outlined here for our
specific problem. The Hilbert space <math>\mathcal{H}</math>\ is decomposed into three
subspaces. The incoming space, denoted by <math>D_{-},</math> consists of all waves
travelling towards the plate, either from the left or the right, as
appropriate. The outgoing subspace, denoted by <math>D_{+},</math> consists of all
waves travelling away from the plate, again either to the left or right, as
appropriate. What remains is the scattering space, denoted by <math>K,\ </math>
consisting of the potential and displacement under the plate.
To apply the Lax-Phillips scattering the following conditions are required: <math>
D_{-}<math> and </math>D_{+}</math> must be orthogonal; the incoming subspace must span the
entire space under temporal evolution. For our system, the first condition
follows from the inner product and the second condition follows from the
simple structure of the eigenfunctions of the operator <math>\mathcal{P}.</math> From
the Lax-Phillips scattering theory, since these conditions hold, the
equation of motion for the plate in the absence of incoming waves can be
written
<center><math>
\frac{1}{i}\partial _{t}U=\mathcal{B}U,  (18)
</math></center>
where <math>\mathcal{B}</math> is a non-self-adjoint operator. <math>\mathcal{B}</math> is related
to <math>\mathcal{P}</math> by
<center><math>
e^{i\mathcal{B}t}=P_{K}\left. e^{i\mathcal{P}t}\right| _{K}
</math></center>
where <math>P_{K}</math> is the projection onto the subspace <math>K</math> and <math>\left. {}\right|
_{K}<math> denotes a restriction to </math>K.<math> Therefore </math>e^{i\mathcal{B}t}</math> is the
restricted to <math>K</math> of the evolution of an initial condition which is zero
outside <math>K.</math> It should be noted that the equality in equation (\ref
{nonselffirst}) is in general only true asymptotically. However the
numerical results show for our case we have equality for all times.
From the Lax-Phillips scattering theory, equation (18) can
be solved by finding the eigenvalues (sometimes referred to as scattering
frequencies or resonances) and eigenfunctions of <math>\mathcal{B}</math>. The
eigenvalues of <math>\mathcal{B}</math> occur at the singularities of the analytic
extension to <math>\mathbb{C}</math> of the scattering matrix, <math>\mathbf{S}(\lambda ).</math>
This is given by
<center><math>
\mathbf{S}(\lambda )=\left(
\begin{matrix}{cc}
S_{11}\left( \lambda \right) & S_{12}\left( \lambda \right) \\
S_{21}\left( \lambda \right) & S_{22}\left( \lambda \right)
\end{matrix}
\right)  (19)
</math></center>
where <math>S_{11}\left( \lambda \right) </math>, <math>S_{12}\left( \lambda \right) ,</math> <math>
S_{21}\left( \lambda \right) ,<math> and </math>S_{22}\left( \lambda \right) </math> are the
scattered wave coefficients found from the single frequency solutions in
section 12. As a consequence of the Lax-Phillips scattering
theory the scattering matrix is unitary for real <math>\lambda </math> and the
singularities must all lie in the upper complex plane (<math>{Im}\left(
\lambda \right) >0).</math> Once the singularities have been found, the
eigenfunctions can be calculated. They are not orthogonal, since <math>\mathcal{B}
</math> is a non-self-adjoint operator, but a biorthogonal system can be formed
using the eigenfunctions of the adjoint operator, <math>\mathcal{B}^{\ast }.</math>
The eigenfunctions of <math>\mathcal{B}</math> are the modes of vibration for the
plate-water system. These modes have a decay as well as an oscillation due
to the radiation of energy into the surrounding water. The frequency of the
oscillation is determined by the real part of the eigenvalue and the rate of
decay is determined by the imaginary part of the eigenvalue.
While the eigenvalues of <math>\mathcal{B}</math> occur precisely at the singularities
of the solution found by a Laplace transformation in time the Lax-Phillips
scattering theory solution has three major advantages over the Laplace
transform solution: the eigenvalues (singularities) can be found using the
scattering matrix; the difficult equations in the Laplace space involving
the initial condition do not need to be solved; the contribution of the
singularity (the residue) can be found directly from the inner product of
the initial condition with the corresponding eigenfunction of the adjoint
operator, <math>\mathcal{B}^{\ast }.</math>
===Finding the Singularities of the Scattering Matrix===
While the analytic extension of the scattering matrix is straightforward,
the linear system (15) is solved for complex <math>\lambda ,</math>
finding the singularities of the scattering matrix is non-trivial<math>.</math> The
difficulty lies in the fact that we must search the complex plane for the
singularities with no ''a priori''  knowledge about their location. We use
a complex search algorithm based contour integration. The determinant of the
scattering matrix is integrated around the contour of a region of the
complex plane. If the value of this integral is zero, then the region is
assumed to contain no singularities and the search is terminated (the
possibility that the contribution of two or more singularities might cancel
can be treated by considering further integrals, such as the variation of
the argument around the contour). If the value of the integral is not zero,
then the region must contain singularities and it is then divided into
subregions and the search is repeated. Once the singularities have been
located sufficiently well they are used as seeds for Newton's method and
found to high accuracy.
If the eigenvalues have to be found for different parameter values then a
homotopy, or continuation, method can be used, which avoids the slow complex
search method. This method uses the known locations of the eigenvalues for
one parameter value to determine the eigenvalues for a new parameter value
by taking sufficiently small steps that Newton's method can be used with the
previous solution as a seed. Unfortunately, a homotopy method requires the
solution of the eigenvalues for at least one parameter value as an initial
seed and this must be accomplished by a complex search algorithm.
The position of the eigenvalues for <math>\beta =2\times 10^{4}</math> and <math>b=50</math> are
shown in Figure \ref{spectrum}. They are denoted by <math>\lambda _{n},</math> where <math>
n\in \mathbb{Z}<math>, and ordered by increasing real part, with </math>n=0</math>
corresponding the eigenvalue with smallest absolute real part. From the
picture and on physical grounds, it seems likely that there exist
asymptotics for the eigenvalues, however this theory is not developed here.
===Eigenfunctions===
The eigenfunctions of <math>\mathcal{B}</math> associated with the eigenvalue <math>\lambda
_{n}<math> are denoted by </math>\Phi ^{+}(\lambda _{n},x),<math> and those of </math>\mathcal{B}
^{\ast }<math> (the adjoint of </math>\mathcal{B})<math> associated with the eigenvalue </math>
\lambda _{n}^{\ast }<math> are denoted by </math>\hat{\Phi}^{+}\left( \lambda
_{n}^{\ast },x\right) </math>. That is,
<center><math>
\mathcal{B}\Phi ^{+}\left( \lambda _{n},x\right) =\lambda _{n}\Phi
^{+}\left( \lambda _{n},x\right)
</math></center>
and
<center><math>
\mathcal{B}^{\ast }\hat{\Phi}^{+}\left( \lambda _{n}^{\ast },x\right)
=\lambda _{n}^{\ast }\hat{\Phi}^{+}\left( \lambda _{n}^{\ast },x\right) .
</math></center>
The eigenfunction <math>\Phi ^{+}\left( \lambda _{n},x\right) </math> can be written
<center><math>
\Phi ^{+}\left( \lambda _{n},x\right) =\left(
\begin{matrix}{c}
\phi ^{+}\left( \lambda _{n},x\right) \\
i\zeta ^{+}\left( \lambda _{n},x\right)
\end{matrix}
\right) =\left(
\begin{matrix}{c}
\sum\limits_{j=1}^{6}\alpha _{j}e^{\mu _{j}\left( \lambda _{n}\right) x} \\
\sum\limits_{j=1}^{6}-\frac{\alpha _{j}\mu _{j}\left( \lambda _{n}\right)
^{2}}{\lambda _{n}}e^{\mu _{j}\left( \lambda _{n}\right) x}
\end{matrix}
\right)
</math></center>
where <math>\mu _{j}\left( \lambda \right) </math> are the six roots of equation (\ref
{cubicinlambda}) and the coefficients <math>\alpha _{j}</math> are found from the
boundary conditions at the end of the plate (4) and the
following condition. Since the scattering matrix is singular, there are no
incoming wave from either direction. We use the condition that there is no
incoming wave at <math>x=-b</math> and that the outgoing wave is of unit amplitude,
i.e.
<center><math>
\phi ^{+}\left( \lambda _{n},-b\right) =e^{i\lambda _{n}b},\;=\ = \left.
\partial _{x}\phi ^{+}\left( \lambda _{n},x\right) \right| _{x=-b}=i\lambda
_{n}e^{i\lambda _{n}b}.
</math></center>
We do not use the condition that there is no outgoing wave at <math>x=b</math> because
the system will become over determined. Therefore, the coefficients <math>\alpha
_{j}</math> satisfy the linear equation
<center><math>
\mathbf{M}\vec{a}\mathbf{=}\vec{b},
</math></center>
where
<center><math>
\mathbf{M}=\left(
\begin{matrix}{cccccc}
\mu _{1}^{4}e^{-\mu _{1}b} & \mu _{2}^{4}e^{-\mu _{2}b} & \mu
_{3}^{4}e^{-\mu _{3}b} & \mu _{4}^{4}e^{-\mu _{4}b} & \mu _{5}^{4}e^{-\mu
_{5}b} & \mu _{6}^{4}e^{-\mu _{6}b} \\
\mu _{1}^{5}e^{-\mu _{1}b} & \mu _{2}^{5}e^{-\mu _{2}b} & \mu
_{3}^{5}e^{-\mu _{3}b} & \mu _{4}^{5}e^{-\mu _{4}b} & \mu _{5}^{5}e^{-\mu
_{5}b} & \mu _{6}^{5}e^{-\mu _{6}b} \\
\mu _{1}^{4}e^{\mu _{1}b} & \mu _{2}^{4}e^{\mu _{2}b} & \mu _{3}^{4}e^{\mu
_{3}b} & \mu _{4}^{4}e^{\mu _{4}b} & \mu _{5}^{4}e^{\mu _{5}b} & \mu
_{6}^{4}e^{\mu _{6}b} \\
\mu _{1}^{5}e^{\mu _{1}b} & \mu _{2}^{5}e^{\mu _{2}b} & \mu _{3}^{5}e^{\mu
_{3}b} & \mu _{4}^{5}e^{\mu _{4}b} & \mu _{5}^{5}e^{\mu _{5}b} & \mu
_{6}^{5}e^{\mu _{6}b} \\
e^{-\mu _{1}b} & e^{-\mu _{2}b} & e^{-\mu _{3}b} & e^{-\mu _{4}b} & e^{-\mu
_{5}b} & e^{-\mu _{6}b} \\
\mu _{1}e^{-\mu _{1}b} & \mu _{2}e^{-\mu _{2}b} & \mu _{3}e^{-\mu _{3}b} &
\mu _{4}e^{-\mu _{4}b} & \mu _{5}e^{-\mu _{5}b} & \mu _{6}e^{-\mu _{6}b}
\end{matrix}
\right)
</math></center>
and <math>\vec{a}</math> and <math>\vec{b}</math> are given by
<center><math>
\vec{a}=\left(
\begin{matrix}{c}
\alpha _{1} \\
\alpha _{2} \\
\alpha _{3} \\
\alpha _{4} \\
\alpha _{5} \\
\alpha _{6}
\end{matrix}
\right) ,\;\;\;\;\;\vec{b}=\left(
\begin{matrix}{c}
0 \\
0 \\
0 \\
0 \\
e^{i\lambda b} \\
i\lambda e^{i\lambda b}
\end{matrix}
\right) .
</math></center>
The eigenfunctions for the adjoint operator are found similarly.
Figure \ref{eigfunctions} shows the real and imaginary parts of the
eigenfunctions of <math>\mathcal{B}</math> for <math>n=1,</math> <math>3</math>, 5, and <math>7,</math> again with <math>
\beta =2\times 10^{4}<math> and </math>b=50.</math> While the eigenfunctions do not have a
simple shape, increasing oscillation is apparent as <math>n</math> increases.
===Inner products===
A biorthogonal system with respect to the energy inner product (\ref
{energyinnerprod}) is formed by the eigenfunctions of <math>\mathcal{B}</math>, <math>\Phi
^{+}\left( \lambda _{n},x\right) ,<math> and the eigenfunctions of </math>\mathcal{B}
^{\ast },<math> </math>\hat{\Phi}^{+}\left( \lambda _{n},x\right) </math>. To normalise the
biorthogonal system, the inner product of <math>\Phi ^{+}\left( \lambda
_{n},x\right) <math> and </math>\hat{\Phi}^{+}\left( \lambda _{n},x\right) </math> has to be
determined. From the definition of the energy inner product (\ref
{energyinnerprod})
<center><math>\begin{matrix}
\left\langle \Phi \left( \lambda _{m},x\right) ,\hat{\Phi}\left( \lambda
_{n},x\right) \right\rangle _{\mathcal{H}} &=&\int_{-b}^{b}\partial _{x}\phi
^{+}\left( \lambda _{m},x\right) \left( \partial _{x}\hat{\phi}^{+}\left(
\lambda _{n}^{\ast },x\right) \right) ^{\ast }dx  (20) \\
&&+\int_{-b}^{b}(1+P)i\zeta ^{+}\left( \lambda _{m}^{\ast },x\right) \left( i
\hat{\zeta}^{+}\left( \lambda _{n}^{\ast },x\right) \right) ^{\ast }dx.
\nonumber
\end{matrix}</math></center>
The two integrals in (20) are considered separately. The
first is
<center><math>\begin{matrix}
&&\int_{-b}^{b}\partial _{x}\phi ^{+}\left( \lambda _{m},x\right) \left(
\partial _{x}\hat{\phi}^{+}\left( \lambda _{n}^{\ast },x\right) \right)
^{\ast }dx \\
&=&\int_{-b}^{b}\left( \sum_{j=1}^{6}\mu _{j}\left( \lambda _{m}\right)
\alpha _{j}e^{\mu _{j}\left( \lambda _{m}\right) x}\right) \left(
\sum_{k=1}^{6}\mu _{k}\left( \lambda _{m}\right) \alpha _{k}e^{\mu
_{k}\left( \lambda _{n}\right) x}\right) dx \\
&=&\sum_{j=1}^{6}\sum_{k=1}^{6}\int_{-b}^{b}-\mu _{j}\left( \lambda
_{m}\right) ^{2}\alpha _{j}e^{\mu _{j}\left( \lambda _{m}\right) x}\alpha
_{k}e^{\mu _{k}\left( \lambda _{n}\right) x}dx \\
&=&\sum_{j=1}^{6}\sum_{k=1}^{6}-\mu _{j}\left( \lambda _{m}\right)
^{2}\alpha _{j}\alpha _{k}\left( \frac{e^{\left( \mu _{j}\left( \lambda
_{m}\right) +\mu _{k}\left( \lambda _{n}\right) \right) b}-e^{-\left( \mu
_{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) \right) b}
}{\mu _{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) }
\right)  \nonumber
\end{matrix}</math></center>
and the second is
<center><math>\begin{matrix}
&&\int_{-b}^{b}(1+P)i\zeta ^{+}\left( \lambda _{m}^{\ast },x\right) \left( i
\hat{\zeta}^{+}\left( \lambda _{n}^{\ast },x\right) \right) ^{\ast }dx \\
&=&\int_{-b}^{b}\left( \sum_{j=1}^{6}-\frac{\alpha _{j}}{\lambda _{m}}
\left( \mu _{j}\left( \lambda _{m}\right) ^{2}+\beta \mu _{j}\left( \lambda
_{m}\right) ^{6}\right) e^{\mu _{j}\left( \lambda _{m}\right) x}\right)
\left( \sum_{k=1}^{6}-\frac{\alpha _{k}}{\lambda _{n}}\mu _{k}\left(
\lambda _{n}\right) ^{2}e^{\mu _{k}\left( \lambda _{n}\right) x}\right) dx \\
&=&\sum_{j=1}^{6}\sum_{k=1}^{6}\int_{-b}^{b}\frac{\alpha _{j}\alpha _{k}}{
\lambda _{m}\lambda _{n}}\left( \mu _{j}\left( \lambda _{m}\right)
^{2}+\beta \mu _{j}\left( \lambda _{m}\right) ^{6}\right) \mu _{k}\left(
\lambda _{n}\right) ^{2}e^{\mu _{j}\left( \lambda _{m}\right) x}e^{\mu
_{k}\left( \lambda _{n}\right) x}dx \\
&=&\sum_{j=1}^{6}\sum_{k=1}^{6}\frac{\alpha _{j}\alpha _{k}}{\lambda
_{m}\lambda _{n}}\left( \mu _{j}\left( \lambda _{m}\right) ^{2}+\beta \mu
_{j}\left( \lambda _{m}\right) ^{6}\right) \mu _{k}\left( \lambda
_{n}\right) ^{2}\left( \frac{e^{\left( \mu _{j}\left( \lambda _{m}\right)
+\mu _{k}\left( \lambda _{n}\right) \right) b}-e^{-\left( \mu _{j}\left(
\lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) \right) b}}{\mu
_{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) }\right)
\\
&=&\sum_{j=1}^{6}\sum_{k=1}^{6}\frac{\alpha _{j}\alpha _{k}}{\lambda
_{m}\lambda _{n}}\left( -\lambda _{m}^{2}\right) \mu _{k}\left( \lambda
_{n}\right) ^{2}\left( \frac{e^{\left( \mu _{j}\left( \lambda _{m}\right)
+\mu _{k}\left( \lambda _{n}\right) \right) b}-e^{-\left( \mu _{j}\left(
\lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) \right) b}}{\mu
_{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) }\right) .
\nonumber
\end{matrix}</math></center>
Therefore the calculation of the inner product in equation (\ref
{innerprodall}) does not require numerical integration.
===Solution===
By solving (18) using the eigenfunctions of <math>\mathcal{B}</math>
and <math>\mathcal{B}^{\ast }</math> we find the evolution of the plate, from an
initial displacement <math>U_{0}(x),</math> is
<center><math>
U\left( x,t\right) =\sum_{n=-\infty }^{\infty }e^{i\lambda _{n}t}\frac{
\left\langle U_{0}\left( x\right) ,\hat{\Phi}\left( \lambda _{n},x\right)
\right\rangle _{\mathcal{H}}}{\left\langle \Phi \left( \lambda _{n},x\right)
,\hat{\Phi}\left( \lambda _{n},x\right) \right\rangle _{\mathcal{H}}}\Phi
\left( \lambda _{n},x\right) .  (21)
</math></center>
The inner product of <math>U_{0}</math> with the eigenfunction <math>\hat{\Phi}\left(
\lambda _{n},x\right) </math> is the only quantity left to compute in (\ref
{evolutionB}). This inner product is written
<center><math>
\left\langle U_{0}\left( x\right) ,\hat{\Phi}\left( \lambda _{n},x\right)
\right\rangle _{\mathcal{H}}=\sum\limits_{j=1}^{6}\left(
\begin{matrix}{c}
\alpha _{j}\dint_{-b}^{b}-\mu _{j}\left( \lambda _{n}\right) ^{2}e^{\mu
_{j}\left( \lambda _{n}\right) x}\phi _{0}\left( x\right) dx+\left. \alpha
_{j}\mu _{j}\left( \lambda _{n}\right) e^{\mu _{j}\left( \lambda _{n}\right)
x}\phi _{0}\left( x\right) \right| _{-b}^{b} \\
\alpha _{j}\lambda _{n}\dint_{-b}^{b}e^{\mu _{j}\left( \lambda _{n}\right)
x}i\zeta _{0}\left( x\right) dx
\end{matrix}
\right)  (22)
</math></center>
and the integrals must be evaluated by numerical integration. The solutions
calculated using the Lax-Phillips scattering theory are identical to those
found using the self-adjoint operator method and for this reason no further
figures are shown.\pagebreak
==Conclusions==
The spectral theory of a linear thin plate floating on shallow water has
been derived. Two spectral-theory solutions have been presented which
determine the time-dependent motion of the thin plate. The first method was
based on self-adjoint operator theory and the second method was based on the
Lax-Phillips scattering. The self-adjoint method solved both the wave
forcing and the free plate problem while the Lax-Phillips method only solved
for a free plate. The eigenfunctions for the self-adjoint method are
orthogonal and the eigenvalues are continuous and consist of all <math>\mathbb{R}
, </math> which makes the calculation of the eigenvalues trivial. The Lax-Phillips
method has discrete eigenvalues which must be calculated numerically and the
system of eigenfunctions is biorthogonal. The advantage of the Lax-Phillips
method is that the modes of vibration of the plate-water system and their
frequency and rate of decay are found. While the relative speeds of the two
methods depends of the exact way in which they are implemented, the
Lax-Phillips method should be considerably faster if the eigenvalues have
been determined.
The development of a spectral theory for more complicated hydroelastic
problems remains a major challenge. While this theory must be more
complicated than that presented here, many features can be expected to
remain. For example, [[ohmatsuVLFS]] has shown that the single frequency
solutions can be used to solve certain time-dependent problems and \cite
{Hazard} have shown that modes, in which the solution can be expanded, exist
for other hydroelastic systems.
{\Large Acknowledgments}
\begin{acknowledgment}
I would like to thank the anonymous reviewers, Dr. Kathy Ruggerio, and Prof.
James Sneyd for their very helpful comments. Also, Prof. Boris Pavlov for
explaining the Lax-Phillips scattering. \pagebreak
\end{acknowledgment}
\bibliographystyle{jfm}
\bibliography{mike,others}
\pagebreak
\begin{center}
{\huge Figure Captions}
\end{center}
\textsc{Figure} 1. Schematic diagram of a thin plate floating on shallow
water and the coordinates and dimensions of the problem.
\textsc{Figure} 2. The evolution of the displacement due to a pulse
travelling to the right for the times shown. The plate occupies the region <math>
-50\leq x\leq 50<math> and is shown by the bold line. </math>\beta =2\times
10^{4},b=50. </math>
\textsc{Figure} 3. The evolution of the potential due to a pulse travelling
to the right for the times shown. The plate occupies the region <math>-50\leq
x\leq 50<math> and is shown by the bold line. </math>\beta =2\times 10^{4},b=50.</math>
\textsc{Figure} 4. The evolution of the displacement for a plate released at
<math>t=0</math> for the times shown. The plate occupies the region <math>-50\leq x\leq 50</math>
and is shown by the bold line. <math>\beta =2\times 10^{4},b=50.</math>
\textsc{Figure} 5. The evolution of the displacement for a plate released at
<math>t=0</math> for the times shown. The plate occupies the region <math>-50\leq x\leq 50</math>
and is shown by the bold line. <math>\beta =2\times 10^{4},b=50.</math>
\textsc{Figure} 6. The location of the first 19 eigenvalues <math>\lambda _{n}</math>
of <math>\mathcal{B}</math> for <math>\beta =2\times 10^{4},\;b=50.</math>
\textsc{Figure} 7. The real (solid) and imaginary (dashed) parts of the
resonance eigenfunctions for <math>n=1,</math> <math>3</math>, 5, and <math>7</math> as shown. <math>\beta
=2\times 10^{4},\;b=50.</math>
\pagebreak
\FRAME{ftFU}{5.5434in}{1.6544in}{0pt}{\Qcb{}}{\Qlb{shallow}}{shallow.eps}{
\special{language "Scientific Word";type "GRAPHIC";display "ICON";valid_file
"F";width 5.5434in;height 1.6544in;depth 0pt;original-width
6.9297in;original-height 1.6544in;cropleft "0";croptop "1";cropright
"1";cropbottom "0";filename 'shallow.eps';file-properties "XNPEU";}}
Some nonsense\pagebreak
\FRAME{ftF}{4.6164in}{3.6296in}{0pt}{}{\Qlb{incomingspecplot1}}{
incomingspecplot1.eps}{\special{language "Scientific Word";type
"GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width
4.6164in;height 3.6296in;depth 0pt;original-width 6.8632in;original-height
5.3869in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename
'incomingspecplot1.eps';file-properties "XNPEU";}}
\FRAME{ftF}{4.5602in}{3.6288in}{0pt}{}{\Qlb{incomingspecplot2}}{
incomingspecplot2.eps}{\special{language "Scientific Word";type
"GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width
4.5602in;height 3.6288in;depth 0pt;original-width 6.7793in;original-height
5.386in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename
'incomingspecplot2.eps';file-properties "XNPEU";}}
\FRAME{fthFU}{4.6155in}{3.6357in}{0pt}{\Qcb{{}}}{\Qlb{spectral1}}{
spectral1.eps}{\special{language "Scientific Word";type
"GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width
4.6155in;height 3.6357in;depth 0pt;original-width 6.8493in;original-height
5.3852in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename
'spectral1.eps';file-properties "XNPEU";}}
\FRAME{fthFU}{4.5887in}{3.608in}{0pt}{\Qcb{{}}}{\Qlb{spectral2}}{
spectral2.eps}{\special{language "Scientific Word";type
"GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width
4.5887in;height 3.608in;depth 0pt;original-width 6.8493in;original-height
5.3852in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename
'spectral2.eps';file-properties "XNPEU";}}
\FRAME{ftFU}{4.5982in}{3.7853in}{0pt}{\Qcb{{}}}{\Qlb{spectrum}}{spectrum.eps
}{\special{language "Scientific Word";type "GRAPHIC";maintain-aspect-ratio
TRUE;display "ICON";valid_file "F";width 4.5982in;height 3.7853in;depth
0pt;original-width 6.8165in;original-height 5.5486in;cropleft "0";croptop
"1";cropright "1";cropbottom "0";filename 'spectrum.eps';file-properties
"XNPEU";}}
\FRAME{ftbpFU}{4.6164in}{3.7377in}{0pt}{\Qcb{{}}}{\Qlb{eigfunctions}}{
eigfunctions.eps}{\special{language "Scientific Word";type
"GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width
4.6164in;height 3.7377in;depth 0pt;original-width 6.845in;original-height
5.5365in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename
'eigfunctions.eps';file-properties "XNPEU";}}

Revision as of 07:48, 25 July 2006

Michael H Meylan, Spectral Solution of Time Dependent Shallow Water Hydroelasticity, J. of Fluid Mechanics, 454, pp 387-402, 2002.

Time-Dependent Linear Water Wave problem of a Floating Elastic Plate on Shallow Depth. The solution was found using a Generalised Eigenfunction Expansion and as a sum over Scattering Frequencies.


Spectral Solution of Time Dependent Shallow Water Hydroelasticity

\author{M\ls I\ls C\ls H\ls A\ls E\ls L\ns H.\ns M\ls E\ls Y\ls L\ls A\ls N[math]\displaystyle{ ^{1} }[/math]} ?? and in revised form ??


The spectral theory of a thin plate floating on shallow water is derived and used to solve the time-dependent motion. This theory is based on an energy inner product in which the evolution operator becomes unitary. Two solution methods are presented. In the first, the solution is expanded in the eigenfunctions of a self-adjoint operator, which are the incoming wave solutions for a single frequency. In the second, the scattering theory of Lax-Phillips is used. The Lax-Phillips scattering solution is suitable for calculating only the free motion of the plate. However, it determines the modes of vibration of the plate-water system. These modes, both oscillate and decay, are found by a complex search algorithm based contour integration. As well as an application to modelling floating runways, the spectral-theory for a floating thin plate on shallow water is a solvable model for more complicated hydroelastic systems.


Introduction

Hydroelasticity is the study of immersed or floating elastic bodies in a fluid. It has a wide range of applications including very large floating structures (Kashiwagihydro98), ships (Bishop), breakwaters ( Stoker) and sea ice (Squire_Review). One of the best studied hydroelastic models is the linear floating thin plate (Newman_deform, ohmatsuVLFS, Kagemoto97, and Kashiwagihydro98) because it models many physical systems, such as a floating runway or an ice floe ( jgrfloecirc).

The time-dependence in linear hydroelastic problems is usually removed by solving for a single frequency which we will refer to as the \emph{single frequency solution}. The solution is normally found by expanding the elastic body motion in the free modes of vibration and solving the fluid equations using a Green function (Bishop). This is analogous to solving for a rigid floating body using the six rigid modes (Sarp_Isa). While alternative methods have been developed (Kashiwagibspline, \cite {ohmatsuVLFS}, and Kagemoto97), these are based on exploiting some property such as a simple geometry or high relative stiffness.

In contrast to the single frequency solutions, solving time-dependent linear hydroelastic systems remains a major challenge. Kashiwagitime and Endotime have developed a time-stepping procedure; however, this method results in error growth in time. Since the problem is linear it is solvable by a spectral method which eliminates the long-time growth of errors. Furthermore, such a method provides information about the behaviour of the solution, such as the decay constant of the motion. However, the spectral theory for linear hydroelasticity has not been developed. For this reason, spectral type solutions such as Ohmatsutimedep, based on a Fourier expansion of the single frequency solution, have only solved the problem in restricted circumstances.

A floating thin plate on shallow water is considered in this paper. This problem has been chosen for the following reasons: while the single frequency solution is straightforward (Stoker), the time-dependent solution has never been calculated; recently OhkusuNamba, \cite {OhkusuISOPE}, and Ertekinshallow1999 used it to model a floating runway; the spectral-theory developed here is a solvable model for more complex hydroelastic systems.

The spectral theory for a thin plate on shallow water is based on an inner product which gives the energy of the plate-water system. With respect to this inner product the evolution operator becomes unitary. Two different solution methods are derived from this spectral theory. The first method is based on an expansion of the solution in eigenfunctions of a self-adjoint operator. These are the single frequency solutions. The second solution method is based on the scattering theory of Lax-Phillips (laxphilips ). It provides the solution in terms of a countable number of modes which have both an oscillation and a decay. These modes are important as they can be used to characterise the response of the system. With the exception of Hazard, they have not been investigated for hydroelastic systems.

Formulation: A Thin Plate on Shallow Water

Figure \ref{shallow} shows a schematic diagram of the problem. The plate is infinite in the [math]\displaystyle{ y }[/math] direction, so that only the [math]\displaystyle{ x }[/math] and [math]\displaystyle{ z }[/math] directions are considered. The [math]\displaystyle{ x }[/math] direction is horizontal, the positive [math]\displaystyle{ z }[/math] axis points vertically up, and the plate covers the region [math]\displaystyle{ -b\leqslant x\leqslant b. }[/math] The water is of uniform depth [math]\displaystyle{ h }[/math] which is small[math]\displaystyle{ \ }[/math]enough that the water may be approximated as shallow (Stoker). The amplitudes are assumed small enough that the linear theory is appropriate, and the plate is sufficiently thin that the shallow draft approximation may be made (\cite {OhkusuISOPE}). The solution could be extended to waves incident at an angle on a infinite two dimensional plate, as described in OhkusuISOPE, but to keep the treatment straightforward this is not done here.

The mathematical description of the problem follows from Stoker. The kinematic condition is

[math]\displaystyle{ \partial _{t}\zeta =-h\partial _{x}^{2}\phi , (1) }[/math]

where [math]\displaystyle{ \phi }[/math] is the velocity potential of the water (averaged over the depth) and [math]\displaystyle{ \zeta }[/math] is the displacement of the water surface or the plate (from the shallow draft approximation). The equation derived by equating the pressure at the free surface is

[math]\displaystyle{ -\rho g\zeta -\rho \partial _{t}\phi =\left\{ \begin{matrix}{c} 0,\;\;x\notin (-b,b), \\ D\partial _{x}^{4}\zeta +\rho ^{\prime }d\partial _{t}^{2}\zeta ,\;\;x\in (-b,b), \end{matrix} \right. (2) }[/math]

where [math]\displaystyle{ D }[/math] is the bending rigidity of the plate per unit length, [math]\displaystyle{ \rho }[/math] is the density of water, [math]\displaystyle{ \rho ^{\prime } }[/math] is the average density of the plate, [math]\displaystyle{ g }[/math] is the acceleration due to gravity, and [math]\displaystyle{ d }[/math] is the thickness of the plate [math]\displaystyle{ . }[/math] At the ends of the plate the free edge boundary conditions

[math]\displaystyle{ \lim_{x\downarrow -b}\partial _{x}^{2}\zeta =\lim_{x\uparrow b}\partial _{x}^{2}\zeta =\lim_{x\downarrow -b}\partial _{x}^{3}\zeta =\lim_{x\uparrow b}\partial _{x}^{3}\zeta =0 (3) }[/math]

are applied since these are common in offshore engineering applications ( OhkusuISOPE). However the theory which will be developed applies equally to any of the energy-conserving edge conditions such as clamped or pinned and there is no need for the boundary conditions to be symmetric. Equation (3) gives the following implied boundary conditions for [math]\displaystyle{ \phi }[/math]

[math]\displaystyle{ \lim_{x\downarrow -b}\partial _{x}^{4}\phi =\lim_{x\uparrow b}\partial _{x}^{4}\phi =\lim_{x\downarrow -b}\partial _{x}^{5}\phi =\lim_{x\uparrow b}\partial _{x}^{5}\phi =0 (4) }[/math]

which will be used subsequently.

Non-dimensional variables are now introduced. The space variables are non-dimensionalised using the water depth [math]\displaystyle{ h, }[/math] and the time variables are non-dimensionalised using [math]\displaystyle{ \sqrt{h/g} }[/math]. The non-dimensional variables are

[math]\displaystyle{ \bar{x}=\frac{x}{h},\;\bar{t}=t\sqrt{\frac{g}{h}},\;\bar{\zeta}=\frac{\zeta }{h},\;\mathrm{and}\;\bar{\phi}=\frac{\phi }{h^{2}\sqrt{g/h}}. }[/math]

In these new variables, (1) and (2) become

[math]\displaystyle{ \partial _{\bar{t}}\bar{\zeta}=-\partial _{\bar{x}}^{2}\bar{\phi} (5) }[/math]

and

[math]\displaystyle{ -\bar{\zeta}-\partial _{\bar{t}}\bar{\phi}=\left\{ \begin{matrix}{c} 0,\;\;\bar{x}\notin (-\bar{b},\bar{b}), \\ \beta \partial _{\bar{x}}^{4}\bar{\zeta}+\gamma \partial _{\bar{t}}^{2}\bar{ \zeta},\;\;\bar{x}\in (-\bar{b},\bar{b}), \end{matrix} \right. (6) }[/math]

where [math]\displaystyle{ \beta }[/math] and [math]\displaystyle{ \gamma }[/math] are

[math]\displaystyle{ \beta =\frac{D}{\rho gh^{4}}\;\;\mathrm{and\ \ }\gamma =\frac{\rho ^{\prime }d}{\rho h}. }[/math]

For clarity the overbar is dropped from now on.

The main change in extending the formulation to water of finite depth is that the velocity potential will be governed by Laplace's equation. This makes the solution of the problem much more computationally demanding since Laplace's equation must be solved by a numerical method, for example the boundary element method. Furthermore, the extension of the spectral theory, which will be developed here for shallow water, to water of finite depth is non-trivial and remains a subject for further work.

Neglecting the inertia term.

It can be assumed that [math]\displaystyle{ \left| \gamma \partial _{t}^{2}\zeta \right| \ll \left| \zeta \right| }[/math] for the following reasons (OhkusuISOPE). If we consider a mode of the displacement [math]\displaystyle{ \zeta =ae^{i\omega t} }[/math] (where [math]\displaystyle{ a }[/math] is the amplitude) then [math]\displaystyle{ \partial _{t}^{2}\zeta =-\omega ^{2}ae^{i\omega t}. }[/math] For each frequency, [math]\displaystyle{ \omega , }[/math] there is a corresponding wavelength [math]\displaystyle{ \lambda =2\pi /\omega . }[/math] In the non-dimensional variables the wave speed and water depth are both unity. Since the water is shallow the wavelength [math]\displaystyle{ \lambda \gg 1\lt math\gt and thus }[/math]\omega \ll 1.</math> It follows that any shallow water mode must satisfy [math]\displaystyle{ \left| \partial _{t}^{2}\zeta \right| \ll \left| \zeta \right| . }[/math] Also, [math]\displaystyle{ \gamma \ll 1 }[/math] since [math]\displaystyle{ \rho ^{\prime }\lt \rho }[/math] (otherwise the plate would sink) and [math]\displaystyle{ d\ll h }[/math] (otherwise the submergence of the plate would not be negligible). Therefore,

[math]\displaystyle{ \left| \gamma \partial _{t}^{2}\zeta \right| \ll \left| \zeta \right| , }[/math]

and we assume in what follows that the inertia, [math]\displaystyle{ \gamma \partial _{t}^{2}\zeta , }[/math] is zero.

It should be noted that the inclusion of the inertia term in the spectral theory which will be developed is difficult because it introduces a time dependence in the energy inner product.

The Energy Inner Product

While equations (5) and (6) are not complicated they cannot be solved in a simple manner. It is not possible to Fourier transform in space because of the spatial discontinuity of the differential equations. The Weiner-Hopf technique cannot be used because the discontinuities divide the space into three regions. A Laplace transformation in time can be applied but this leads to non-trivial equations involving a spatially discontinuous differential equation subject to arbitrary initial conditions. However, straightforward solutions can be derived using spectral theory.

The spectral-theory solution of equations (5) and (\ref {displacement2}) is based on the spectral theory for a unitary operator (essentially, an operator is unitary if the adjoint is also the inverse). We therefore require an inner product in which the evolution operator is unitary. This inner product, since the system is conservative, is derived from the energy. The potential and displacement both contribute to this energy and we combine them in a two component vector, [math]\displaystyle{ U\left( x,t\right) }[/math], given by

[math]\displaystyle{ U\left( x,t\right) =\left( \begin{matrix}{c} \phi (x,t) \\ i\zeta (x,t) \end{matrix} \right) . (7) }[/math]

The energy consists of the kinetic energy of the water ([math]\displaystyle{ \propto \left| \phi _{t}^{2}\right| \lt math\gt ), the potential energy of the water ( }[/math]\propto \left| \phi ^{2}\right| </math>), and the energy of the plate. The energy inner product for the two vectors

[math]\displaystyle{ U_{1}=\left( \begin{matrix}{c} \phi _{1} \\ i\zeta _{1} \end{matrix} \right) \;\;\mathrm{and\ \ }U_{2}=\left( \begin{matrix}{c} \phi _{2} \\ i\zeta _{2} \end{matrix} \right) }[/math]

is given by

[math]\displaystyle{ \left\langle U_{1},U_{2}\right\rangle _{\mathcal{H}}=\left\langle \partial _{x}\phi _{1},\partial _{x}\phi _{2}\right\rangle +\left\langle \left( 1+\beta \left( H\left( x-b\right) -H\left( x+b\right) \right) \partial _{x}^{4}\right) i\zeta _{1},i\zeta _{2}\right\rangle , (8) }[/math]

where [math]\displaystyle{ H }[/math] is the Heaviside function. The subscript [math]\displaystyle{ \mathcal{H} }[/math] is used to denote the special inner product and the angle brackets without the [math]\displaystyle{ \mathcal{H} }[/math] denote the standard inner product, i.e.

[math]\displaystyle{ \left\langle f\left( x\right) ,g\left( x\right) \right\rangle =\int_{-\infty }^{\infty }f\left( x\right) g^{\ast }\left( x\right) dx. }[/math]

We now write (5) and (6) as

[math]\displaystyle{ \begin{matrix} \frac{1}{i}\partial _{t}U &=&\mathcal{P}U (9) \\ U\left( x,t\right) _{t=0} &=&U_{0}\left( x\right) =\left( \begin{matrix}{c} \phi _{0}(x) \\ i\zeta _{0}(x) \end{matrix} \right) \nonumber \end{matrix} }[/math]

where the operator [math]\displaystyle{ \mathcal{P} }[/math] is

[math]\displaystyle{ \mathcal{P=}\left( \begin{matrix}{cc} 0 & 1+\beta \left( H\left( x-b\right) -H\left( x+b\right) \right) \partial _{x}^{4} \\ -\partial _{x}^{2} & 0 \end{matrix} \right) . }[/math]

[math]\displaystyle{ \mathcal{P} }[/math] is self-adjoint with respect to the inner product (\ref {energyinnerprod}) since [math]\displaystyle{ \mathcal{P} }[/math] satisfies

[math]\displaystyle{ \left\langle \mathcal{P}U_{1},U_{2}\right\rangle _{\mathcal{H}}=\left\langle U_{1},\mathcal{P}U_{2}\right\rangle _{\mathcal{H}} }[/math]

from integration by parts and the boundary conditions at the end of the plate (3). We can express the solution to (\ref {selfadjoint2}) as

[math]\displaystyle{ U\left( x,t\right) =e^{i\mathcal{P}t}U_{0} (10) }[/math]

where [math]\displaystyle{ e^{i\mathcal{P}t} }[/math] is a unitary evolution operator.

The self-adjoint solution method(11)

In this section, a solution for the time dependent motion of the plate-water system is developed using the theory of self-adjoint operators. To evaluate equation (10) we require a method to calculate the evolution operator [math]\displaystyle{ e^{i\mathcal{P}t} }[/math]. This can be accomplished by using the eigenfunctions of the operator [math]\displaystyle{ \mathcal{P}, }[/math] which are the single frequency solutions.

Finding the eigenfunctions(12)

Since [math]\displaystyle{ \mathcal{P} }[/math] is self-adjoint, the eigenvalues, [math]\displaystyle{ \lambda , }[/math] must be real and therefore[math]\displaystyle{ \ }[/math]the eigenfunctions of [math]\displaystyle{ \mathcal{P} }[/math] are oscillatory exponentials outside the region of water covered by the plate. Furthermore, since the plate is finite, the spectrum (set of eigenvalues) is the entire real numbers. As is expected for two-component systems, there are two eigenfunctions associated with each eigenvalue [math]\displaystyle{ \lambda }[/math]. We choose incoming waves from the left ([math]\displaystyle{ \Phi ^{\gt }) }[/math] and the right ([math]\displaystyle{ \Phi ^{\lt }) }[/math] of unit amplitude as a basis for the eigenspace since they are the standard single frequency solutions. They have the following asymptotics,

[math]\displaystyle{ \lim_{x\rightarrow -\infty }\Phi ^{\gt }=\left( \begin{matrix}{c} e^{i\lambda x} \\ \lambda e^{i\lambda x} \end{matrix} \right) +S_{12}\left( \begin{matrix}{c} e^{-i\lambda x} \\ \lambda e^{-i\lambda x} \end{matrix} \right) \;\;\mathrm{and\ \ }\lim_{x\rightarrow \infty }\Phi ^{\gt }=S_{11}\left( \begin{matrix}{c} e^{i\lambda x} \\ \lambda e^{i\lambda x} \end{matrix} \right) }[/math]

and

[math]\displaystyle{ \lim_{t\rightarrow -\infty }\Phi ^{\lt }=S_{22}\left( \begin{matrix}{c} e^{-i\lambda x} \\ \lambda e^{-i\lambda x} \end{matrix} \right) \;\;\mathrm{and\ \ }\lim_{x\rightarrow \infty }\Phi ^{\gt }=\left( \begin{matrix}{c} e^{-i\lambda x} \\ \lambda e^{-i\lambda x} \end{matrix} \right) +S_{21}\left( \begin{matrix}{c} e^{i\lambda x} \\ \lambda e^{i\lambda x} \end{matrix} \right) , }[/math]

where [math]\displaystyle{ S_{11} }[/math], [math]\displaystyle{ S_{12,} }[/math] [math]\displaystyle{ S_{21}, }[/math] and [math]\displaystyle{ S_{22} }[/math] are the reflection and transmission coefficients (which must be determined). These eigenfunctions, which are analogous to the Jost solutions of Schrodinger's equation (\cite {Chadan89}), will be used to calculated the time-dependent solution.

We find the eigenfunction [math]\displaystyle{ \Phi ^{\gt }(\lambda ,x) }[/math] by solving (\ref {kinematic2}) and (6) in each region. The two components, [math]\displaystyle{ \phi ^{\gt }\left( \lambda ,x\right) \lt math\gt and }[/math]i\zeta ^{>}\left( \lambda ,x\right) ,</math> are given by

[math]\displaystyle{ \phi ^{\gt }\left( \lambda ,x\right) =\left\{ \begin{matrix}{c} e^{-i\lambda x}+S_{11}\left( \lambda \right) e^{i\lambda x},\;\;x\lt -b, \\ \sum\limits_{j=1}^{6}\alpha _{j}e^{\mu _{j}\left( \lambda \right) x},\;\;-b\lt x\lt b, \\ S_{12}\left( \lambda \right) e^{-i\lambda x},\;\;x\gt b, \end{matrix} \right. (13) }[/math]

and

[math]\displaystyle{ i\zeta ^{\gt }\left( \lambda ,x\right) =\left\{ \begin{matrix}{c} \lambda e^{-i\lambda x}+\lambda S_{11}\left( \lambda \right) e^{i\lambda x},\;\;x\lt -b, \\ \frac{-1}{\lambda }\sum\limits_{j=1}^{6}\mu _{j}\left( \lambda \right) ^{2}\alpha _{j}e^{\mu _{j}\left( \lambda \right) x},\;\;-b\lt x\lt b, \\ \lambda S_{12}\left( \lambda \right) e^{-i\lambda x},\;\;x\gt b, \end{matrix} \right. }[/math]

where the coefficients [math]\displaystyle{ \mu _{j}\left( \lambda \right) }[/math] are the six roots of the equation

[math]\displaystyle{ \beta \mu ^{6}+\mu ^{2}+\lambda ^{2}=0. (14) }[/math]

The values of [math]\displaystyle{ S_{11}\left( \lambda \right) , }[/math] [math]\displaystyle{ S_{12}\left( \lambda \right) ,\lt math\gt and }[/math]\alpha _{j}</math> are found from the boundary conditions (4) and the continuity of [math]\displaystyle{ \phi }[/math] and [math]\displaystyle{ \partial _{x}\phi }[/math] at [math]\displaystyle{ x=\pm b. }[/math] Therefore, to find the eigenfunction [math]\displaystyle{ \Phi ^{\gt }\left( \lambda ,x\right) , }[/math] we solve the 8 by 8 linear system

[math]\displaystyle{ \mathbf{M}\vec{a}\mathbf{=}\vec{b}, (15) }[/math]

where [math]\displaystyle{ \mathbf{M} }[/math] is the matrix

[math]\displaystyle{ \mathbf{M}=\left( \begin{matrix}{cccccccc} \mu _{1}^{4}e^{-\mu _{1}b} & \mu _{2}^{4}e^{-\mu _{2}b} & \mu _{3}^{4}e^{-\mu _{3}b} & \mu _{4}^{4}e^{-\mu _{4}b} & \mu _{5}^{4}e^{-\mu _{5}b} & \mu _{6}^{4}e^{-\mu _{6}b} & 0 & 0 \\ \mu _{1}^{5}e^{-\mu _{1}b} & \mu _{2}^{5}e^{-\mu _{2}b} & \mu _{3}^{5}e^{-\mu _{3}b} & \mu _{4}^{5}e^{-\mu _{4}b} & \mu _{5}^{5}e^{-\mu _{5}b} & \mu _{6}^{5}e^{-\mu _{6}b} & 0 & 0 \\ \mu _{1}^{4}e^{\mu _{1}b} & \mu _{2}^{4}e^{\mu _{2}b} & \mu _{3}^{4}e^{\mu _{3}b} & \mu _{4}^{4}e^{\mu _{4}b} & \mu _{5}^{4}e^{\mu _{5}b} & \mu _{6}^{4}e^{\mu _{6}b} & 0 & 0 \\ \mu _{1}^{5}e^{\mu _{1}b} & \mu _{2}^{5}e^{\mu _{2}b} & \mu _{3}^{5}e^{\mu _{3}b} & \mu _{4}^{5}e^{\mu _{4}b} & \mu _{5}^{5}e^{\mu _{5}b} & \mu _{6}^{5}e^{\mu _{6}b} & 0 & 0 \\ e^{-\mu _{1}b} & e^{-\mu _{2}b} & e^{-\mu _{3}b} & e^{-\mu _{4}b} & e^{-\mu _{5}b} & e^{-\mu _{6}b} & -e^{-i\lambda b} & 0 \\ \mu _{1}e^{-\mu _{1}b} & \mu _{2}e^{-\mu _{2}b} & \mu _{3}e^{-\mu _{3}b} & \mu _{4}e^{-\mu _{4}b} & \mu _{5}e^{-\mu _{5}b} & \mu _{6}e^{-\mu _{6}b} & -i\lambda e^{-i\lambda b} & 0 \\ e^{\mu _{1}b} & e^{\mu _{2}b} & e^{\mu _{3}b} & e^{\mu _{4}b} & e^{\mu _{5}b} & e^{\mu _{6}b} & 0 & -e^{-i\lambda b} \\ \mu _{1}e^{\mu _{1}b} & \mu _{2}e^{\mu _{2}b} & \mu _{3}e^{\mu _{3}b} & \mu _{4}e^{\mu _{4}b} & \mu _{5}e^{\mu _{5}b} & \mu _{6}e^{\mu _{6}b} & 0 & i\lambda e^{-i\lambda b} \end{matrix} \right) }[/math]

and [math]\displaystyle{ \vec{a} }[/math] and [math]\displaystyle{ \vec{b} }[/math] are given by

[math]\displaystyle{ \vec{a}=\left( \begin{matrix}{c} \alpha _{1} \\ \alpha _{2} \\ \alpha _{3} \\ \alpha _{4} \\ \alpha _{5} \\ \alpha _{6} \\ S_{11} \\ S_{12} \end{matrix} \right) ,\;\;\;\;\;\vec{b}=\left( \begin{matrix}{c} 0 \\ 0 \\ 0 \\ 0 \\ e^{i\lambda b} \\ -i\lambda e^{i\lambda b} \\ 0 \\ 0 \end{matrix} \right) . }[/math]

Note that the coefficients [math]\displaystyle{ S_{11} }[/math] and [math]\displaystyle{ S_{12} }[/math] are contained in [math]\displaystyle{ \vec{a}. }[/math]

The eigenfunctions for the wave propagating from the right [math]\displaystyle{ \Phi ^{\lt }\left( \lambda ,x\right) \lt math\gt are found similarly. Since }[/math]S_{11}</math> represents the amplitude of the reflected wave and [math]\displaystyle{ S_{12} }[/math] represents the amplitude of the transmitted wave, conservation of energy requires that [math]\displaystyle{ \left| S_{11}\right| ^{2}+\left| S_{12}\right| ^{2}=1. }[/math] Similarly, since the boundary conditions are symmetric [math]\displaystyle{ S_{22}\left( \lambda \right) =S_{11}\left( \lambda \right) }[/math] and [math]\displaystyle{ S_{12}\left( \lambda \right) =S_{21}\left( \lambda \right) . }[/math]

Solution with the eigenfunctions

Equation (10) can be solved by a generalised Fourier transform based on the eigenfunctions of the operator [math]\displaystyle{ \mathcal{P} }[/math]. The eigenfunctions are orthogonal since [math]\displaystyle{ \mathcal{P} }[/math] is self-adjoint, but they must be normalised. This is accomplished by using the following identity

[math]\displaystyle{ \int_{0}^{\infty }e^{i\left( \lambda _{1}-\lambda _{2}\right) t}dt=\pi \delta \left( \lambda _{1}-\lambda _{2}\right) }[/math]

where [math]\displaystyle{ \delta }[/math] is the Dirac delta function. Therefore

[math]\displaystyle{ \begin{matrix} \left\langle \Phi ^{\gt }\left( x,\lambda _{1}\right) ,\Phi ^{\gt }\left( x,\lambda _{2}\right) \right\rangle _{\mathcal{H}} &=&\pi \delta \left( \lambda _{1}-\lambda _{2}\right) \lambda _{1}^{2}\left( 1+\left| S_{11}\right| ^{2}+\left| S_{12}\right| ^{2}\right) \\ &&+\pi \delta \left( \lambda _{1}-\lambda _{2}\right) \lambda _{1}\lambda _{2}\left( 1+\left| S_{11}\right| ^{2}+\left| S_{12}\right| ^{2}\right) \\ &=&4\pi \delta \left( \lambda _{1}-\lambda _{2}\right) \lambda _{1}^{2}, \nonumber \end{matrix} }[/math]

using the condition [math]\displaystyle{ \left| S_{11}\right| ^{2}+\left| S_{12}\right| ^{2}=1. }[/math] Similarly

[math]\displaystyle{ \left\langle \Phi ^{\lt }\left( x,\lambda _{1}\right) ,\Phi ^{\lt }\left( x,\lambda _{2}\right) \right\rangle _{\mathcal{H}}=4\pi \delta \left( \lambda _{1}-\lambda _{2}\right) \lambda _{1}^{2} }[/math]

and

[math]\displaystyle{ \begin{matrix} \left\langle \Phi ^{\gt }\left( x,\lambda _{1}\right) ,\Phi ^{\lt }\left( x,\lambda _{2}\right) \right\rangle _{\mathcal{H}} &=&2\pi \delta \left( \lambda _{1}-\lambda _{2}\right) \lambda _{1}^{2}\left( S_{11}S_{21}^{\ast }+S_{12}S_{22}^{\ast }\right) \\ &=&0. \nonumber \end{matrix} }[/math]

The generalised Fourier transform which solves the evolution equation (\ref {unitary}) is

[math]\displaystyle{ \begin{matrix} U\left( x,t\right) &=&\int_{-\infty }^{\infty }\left\langle U_{0}\left( x\right) ,\frac{\Phi ^{\gt }\left( x,\lambda \right) }{4\pi \lambda ^{2}} \right\rangle _{\mathcal{H}}\Phi ^{\gt }\left( x,\lambda \right) e^{i\lambda t}d\lambda (16) \\ &&+\int_{-\infty }^{\infty }\left\langle U_{0}\left( x\right) ,\frac{\Phi ^{\lt }\left( x,\lambda \right) }{4\pi \lambda ^{2}}\right\rangle _{\mathcal{H} }\Phi ^{\lt }\left( x,\lambda \right) e^{i\lambda t}d\lambda . \nonumber \end{matrix} }[/math]

Equation (16) is the cornerstone of the approach. The integral in equation (16) can be calculated by the fast Fourier transform while the inner product can calculated by the fast Fourier transform if the initial condition [math]\displaystyle{ U_{0} }[/math] is zero underneath the plate ([math]\displaystyle{ -b\lt x\lt b). }[/math]

Numerical Calculations

The intention of this paper is to develop the solution methods rather than describe the physics of the motion and therefore only a few solutions are presented. From OhkusuISOPE we assume the plate stiffness is [math]\displaystyle{ \beta =2\times 10^{4}\lt math\gt and the plate length is }[/math]b=50</math> throughout. These values are typical for a floating runway. Figures \ref{incomingspecplot1} and \ref {incomingspecplot2} show the displacement and potential, respectively, for a pulse travelling to the right at the times [math]\displaystyle{ t=0, }[/math] 30, 60, 90, 120, 150, 180, 210, and 240. The incoming wave pulse was chosen to be a Gaussian in potential centered at [math]\displaystyle{ x=-125 }[/math] and sufficiently sharp to be negligible under the plate,

[math]\displaystyle{ U_{0}\left( x\right) =\left( \begin{matrix}{c} \phi \left( x\right) \\ i\phi ^{\prime }\left( x\right) \end{matrix} \right) }[/math]

where

[math]\displaystyle{ \phi \left( x\right) =\left\{ \begin{matrix}{c} e^{-\tfrac{(x+125)^{2}}{350}},\;\;x\lt -50, \\ 0,\;\;x\gt -50. \end{matrix} \right. }[/math]

At [math]\displaystyle{ t=0 }[/math] the plate is initially at rest and the wave is to the left of the plate propagating towards it. From [math]\displaystyle{ t=30 }[/math] the wave has reached the plate and the plate begins to undergo a complex bending motion in response to the incoming wave. The response of the plate in turn induces waves in the surrounding water which propagate away from the plate to the left and right. The final picture, [math]\displaystyle{ t=240 }[/math], shows the plate at rest with waves now propagating away from it. The majority of the wave energy has passed under the plate and continues to propagate to the right. However, the shape of the outgoing wave profile is markedly different from the incoming wave profile. Also, there is a significant reflected wave propagating away from the plate to the left.

Figures \ref{spectral1} and \ref{spectral2} show the evolution of the plate from an initial displacement in the absence of wave forcing for the times [math]\displaystyle{ t=0, }[/math] 20, 40, 60, 80, 100, 120, 140, and 160. Only the plate displacement is initially non-zero so that

[math]\displaystyle{ U_{0}\left( x\right) =\left( \begin{matrix}{c} 0 \\ i\zeta \left( x\right) \end{matrix} \right) . }[/math]

Figure \ref{spectral1} shows the motion for the symmetric initial plate displacement

[math]\displaystyle{ \zeta \left( x\right) =e^{-\tfrac{x^{2}}{350}}. }[/math]

As the plate evolves the plate vibrates, straightens, and the amplitude decays. The decay is due to the radiation of energy by the waves generated in the surrounding water. A complex wave train is produced by the plate motion and can be seen propagating away from the plate. Figure \ref {spectral2} shows the motion for the non-symmetric initial plate displacement

[math]\displaystyle{ \zeta \left( x\right) =e^{-\tfrac{\left( x-50\right) ^{2}}{350}}. }[/math]

Again as the plate evolves it straightens, vibrates, and decays and induces waves in the surrounding water.

The Lax-Phillips Scattering Solution Method.

In this section, a solution to the time-dependent motion of the plate-water system is developed using the Lax-Phillips scattering theory (\cite {laxphilips}). This solution method will only solve for an initial condition which is zero outside the plate, i.e. [math]\displaystyle{ U_{0}\left( x\right) =0 }[/math] if [math]\displaystyle{ \left| x\right| \gt b }[/math]. However, it calculates the solution by an expansion in a countable number of modes.

Lax-Phillips Scattering(17)

The Lax-Phillips scattering theory will be briefly outlined here for our specific problem. The Hilbert space [math]\displaystyle{ \mathcal{H} }[/math]\ is decomposed into three subspaces. The incoming space, denoted by [math]\displaystyle{ D_{-}, }[/math] consists of all waves travelling towards the plate, either from the left or the right, as appropriate. The outgoing subspace, denoted by [math]\displaystyle{ D_{+}, }[/math] consists of all waves travelling away from the plate, again either to the left or right, as appropriate. What remains is the scattering space, denoted by [math]\displaystyle{ K,\ }[/math] consisting of the potential and displacement under the plate.

To apply the Lax-Phillips scattering the following conditions are required: [math]\displaystyle{ D_{-}\lt math\gt and }[/math]D_{+}</math> must be orthogonal; the incoming subspace must span the entire space under temporal evolution. For our system, the first condition follows from the inner product and the second condition follows from the simple structure of the eigenfunctions of the operator [math]\displaystyle{ \mathcal{P}. }[/math] From the Lax-Phillips scattering theory, since these conditions hold, the equation of motion for the plate in the absence of incoming waves can be written

[math]\displaystyle{ \frac{1}{i}\partial _{t}U=\mathcal{B}U, (18) }[/math]

where [math]\displaystyle{ \mathcal{B} }[/math] is a non-self-adjoint operator. [math]\displaystyle{ \mathcal{B} }[/math] is related to [math]\displaystyle{ \mathcal{P} }[/math] by

[math]\displaystyle{ e^{i\mathcal{B}t}=P_{K}\left. e^{i\mathcal{P}t}\right| _{K} }[/math]

where [math]\displaystyle{ P_{K} }[/math] is the projection onto the subspace [math]\displaystyle{ K }[/math] and [math]\displaystyle{ \left. {}\right| _{K}\lt math\gt denotes a restriction to }[/math]K.[math]\displaystyle{ Therefore }[/math]e^{i\mathcal{B}t}</math> is the restricted to [math]\displaystyle{ K }[/math] of the evolution of an initial condition which is zero outside [math]\displaystyle{ K. }[/math] It should be noted that the equality in equation (\ref {nonselffirst}) is in general only true asymptotically. However the numerical results show for our case we have equality for all times.

From the Lax-Phillips scattering theory, equation (18) can be solved by finding the eigenvalues (sometimes referred to as scattering frequencies or resonances) and eigenfunctions of [math]\displaystyle{ \mathcal{B} }[/math]. The eigenvalues of [math]\displaystyle{ \mathcal{B} }[/math] occur at the singularities of the analytic extension to [math]\displaystyle{ \mathbb{C} }[/math] of the scattering matrix, [math]\displaystyle{ \mathbf{S}(\lambda ). }[/math] This is given by

[math]\displaystyle{ \mathbf{S}(\lambda )=\left( \begin{matrix}{cc} S_{11}\left( \lambda \right) & S_{12}\left( \lambda \right) \\ S_{21}\left( \lambda \right) & S_{22}\left( \lambda \right) \end{matrix} \right) (19) }[/math]

where [math]\displaystyle{ S_{11}\left( \lambda \right) }[/math], [math]\displaystyle{ S_{12}\left( \lambda \right) , }[/math] [math]\displaystyle{ S_{21}\left( \lambda \right) ,\lt math\gt and }[/math]S_{22}\left( \lambda \right) </math> are the scattered wave coefficients found from the single frequency solutions in section 12. As a consequence of the Lax-Phillips scattering theory the scattering matrix is unitary for real [math]\displaystyle{ \lambda }[/math] and the singularities must all lie in the upper complex plane ([math]\displaystyle{ {Im}\left( \lambda \right) \gt 0). }[/math] Once the singularities have been found, the eigenfunctions can be calculated. They are not orthogonal, since [math]\displaystyle{ \mathcal{B} }[/math] is a non-self-adjoint operator, but a biorthogonal system can be formed using the eigenfunctions of the adjoint operator, [math]\displaystyle{ \mathcal{B}^{\ast }. }[/math]

The eigenfunctions of [math]\displaystyle{ \mathcal{B} }[/math] are the modes of vibration for the plate-water system. These modes have a decay as well as an oscillation due to the radiation of energy into the surrounding water. The frequency of the oscillation is determined by the real part of the eigenvalue and the rate of decay is determined by the imaginary part of the eigenvalue.

While the eigenvalues of [math]\displaystyle{ \mathcal{B} }[/math] occur precisely at the singularities of the solution found by a Laplace transformation in time the Lax-Phillips scattering theory solution has three major advantages over the Laplace transform solution: the eigenvalues (singularities) can be found using the scattering matrix; the difficult equations in the Laplace space involving the initial condition do not need to be solved; the contribution of the singularity (the residue) can be found directly from the inner product of the initial condition with the corresponding eigenfunction of the adjoint operator, [math]\displaystyle{ \mathcal{B}^{\ast }. }[/math]

Finding the Singularities of the Scattering Matrix

While the analytic extension of the scattering matrix is straightforward, the linear system (15) is solved for complex [math]\displaystyle{ \lambda , }[/math] finding the singularities of the scattering matrix is non-trivial[math]\displaystyle{ . }[/math] The difficulty lies in the fact that we must search the complex plane for the singularities with no a priori knowledge about their location. We use a complex search algorithm based contour integration. The determinant of the scattering matrix is integrated around the contour of a region of the complex plane. If the value of this integral is zero, then the region is assumed to contain no singularities and the search is terminated (the possibility that the contribution of two or more singularities might cancel can be treated by considering further integrals, such as the variation of the argument around the contour). If the value of the integral is not zero, then the region must contain singularities and it is then divided into subregions and the search is repeated. Once the singularities have been located sufficiently well they are used as seeds for Newton's method and found to high accuracy.

If the eigenvalues have to be found for different parameter values then a homotopy, or continuation, method can be used, which avoids the slow complex search method. This method uses the known locations of the eigenvalues for one parameter value to determine the eigenvalues for a new parameter value by taking sufficiently small steps that Newton's method can be used with the previous solution as a seed. Unfortunately, a homotopy method requires the solution of the eigenvalues for at least one parameter value as an initial seed and this must be accomplished by a complex search algorithm.

The position of the eigenvalues for [math]\displaystyle{ \beta =2\times 10^{4} }[/math] and [math]\displaystyle{ b=50 }[/math] are shown in Figure \ref{spectrum}. They are denoted by [math]\displaystyle{ \lambda _{n}, }[/math] where [math]\displaystyle{ n\in \mathbb{Z}\lt math\gt , and ordered by increasing real part, with }[/math]n=0</math> corresponding the eigenvalue with smallest absolute real part. From the picture and on physical grounds, it seems likely that there exist asymptotics for the eigenvalues, however this theory is not developed here.

Eigenfunctions

The eigenfunctions of [math]\displaystyle{ \mathcal{B} }[/math] associated with the eigenvalue [math]\displaystyle{ \lambda _{n}\lt math\gt are denoted by }[/math]\Phi ^{+}(\lambda _{n},x),[math]\displaystyle{ and those of }[/math]\mathcal{B} ^{\ast }[math]\displaystyle{ (the adjoint of }[/math]\mathcal{B})[math]\displaystyle{ associated with the eigenvalue }[/math] \lambda _{n}^{\ast }[math]\displaystyle{ are denoted by }[/math]\hat{\Phi}^{+}\left( \lambda _{n}^{\ast },x\right) </math>. That is,

[math]\displaystyle{ \mathcal{B}\Phi ^{+}\left( \lambda _{n},x\right) =\lambda _{n}\Phi ^{+}\left( \lambda _{n},x\right) }[/math]

and

[math]\displaystyle{ \mathcal{B}^{\ast }\hat{\Phi}^{+}\left( \lambda _{n}^{\ast },x\right) =\lambda _{n}^{\ast }\hat{\Phi}^{+}\left( \lambda _{n}^{\ast },x\right) . }[/math]

The eigenfunction [math]\displaystyle{ \Phi ^{+}\left( \lambda _{n},x\right) }[/math] can be written

[math]\displaystyle{ \Phi ^{+}\left( \lambda _{n},x\right) =\left( \begin{matrix}{c} \phi ^{+}\left( \lambda _{n},x\right) \\ i\zeta ^{+}\left( \lambda _{n},x\right) \end{matrix} \right) =\left( \begin{matrix}{c} \sum\limits_{j=1}^{6}\alpha _{j}e^{\mu _{j}\left( \lambda _{n}\right) x} \\ \sum\limits_{j=1}^{6}-\frac{\alpha _{j}\mu _{j}\left( \lambda _{n}\right) ^{2}}{\lambda _{n}}e^{\mu _{j}\left( \lambda _{n}\right) x} \end{matrix} \right) }[/math]

where [math]\displaystyle{ \mu _{j}\left( \lambda \right) }[/math] are the six roots of equation (\ref {cubicinlambda}) and the coefficients [math]\displaystyle{ \alpha _{j} }[/math] are found from the boundary conditions at the end of the plate (4) and the following condition. Since the scattering matrix is singular, there are no incoming wave from either direction. We use the condition that there is no incoming wave at [math]\displaystyle{ x=-b }[/math] and that the outgoing wave is of unit amplitude, i.e.

[math]\displaystyle{ \phi ^{+}\left( \lambda _{n},-b\right) =e^{i\lambda _{n}b},\;=\ = \left. \partial _{x}\phi ^{+}\left( \lambda _{n},x\right) \right| _{x=-b}=i\lambda _{n}e^{i\lambda _{n}b}. }[/math]

We do not use the condition that there is no outgoing wave at [math]\displaystyle{ x=b }[/math] because the system will become over determined. Therefore, the coefficients [math]\displaystyle{ \alpha _{j} }[/math] satisfy the linear equation

[math]\displaystyle{ \mathbf{M}\vec{a}\mathbf{=}\vec{b}, }[/math]

where

[math]\displaystyle{ \mathbf{M}=\left( \begin{matrix}{cccccc} \mu _{1}^{4}e^{-\mu _{1}b} & \mu _{2}^{4}e^{-\mu _{2}b} & \mu _{3}^{4}e^{-\mu _{3}b} & \mu _{4}^{4}e^{-\mu _{4}b} & \mu _{5}^{4}e^{-\mu _{5}b} & \mu _{6}^{4}e^{-\mu _{6}b} \\ \mu _{1}^{5}e^{-\mu _{1}b} & \mu _{2}^{5}e^{-\mu _{2}b} & \mu _{3}^{5}e^{-\mu _{3}b} & \mu _{4}^{5}e^{-\mu _{4}b} & \mu _{5}^{5}e^{-\mu _{5}b} & \mu _{6}^{5}e^{-\mu _{6}b} \\ \mu _{1}^{4}e^{\mu _{1}b} & \mu _{2}^{4}e^{\mu _{2}b} & \mu _{3}^{4}e^{\mu _{3}b} & \mu _{4}^{4}e^{\mu _{4}b} & \mu _{5}^{4}e^{\mu _{5}b} & \mu _{6}^{4}e^{\mu _{6}b} \\ \mu _{1}^{5}e^{\mu _{1}b} & \mu _{2}^{5}e^{\mu _{2}b} & \mu _{3}^{5}e^{\mu _{3}b} & \mu _{4}^{5}e^{\mu _{4}b} & \mu _{5}^{5}e^{\mu _{5}b} & \mu _{6}^{5}e^{\mu _{6}b} \\ e^{-\mu _{1}b} & e^{-\mu _{2}b} & e^{-\mu _{3}b} & e^{-\mu _{4}b} & e^{-\mu _{5}b} & e^{-\mu _{6}b} \\ \mu _{1}e^{-\mu _{1}b} & \mu _{2}e^{-\mu _{2}b} & \mu _{3}e^{-\mu _{3}b} & \mu _{4}e^{-\mu _{4}b} & \mu _{5}e^{-\mu _{5}b} & \mu _{6}e^{-\mu _{6}b} \end{matrix} \right) }[/math]

and [math]\displaystyle{ \vec{a} }[/math] and [math]\displaystyle{ \vec{b} }[/math] are given by

[math]\displaystyle{ \vec{a}=\left( \begin{matrix}{c} \alpha _{1} \\ \alpha _{2} \\ \alpha _{3} \\ \alpha _{4} \\ \alpha _{5} \\ \alpha _{6} \end{matrix} \right) ,\;\;\;\;\;\vec{b}=\left( \begin{matrix}{c} 0 \\ 0 \\ 0 \\ 0 \\ e^{i\lambda b} \\ i\lambda e^{i\lambda b} \end{matrix} \right) . }[/math]

The eigenfunctions for the adjoint operator are found similarly.

Figure \ref{eigfunctions} shows the real and imaginary parts of the eigenfunctions of [math]\displaystyle{ \mathcal{B} }[/math] for [math]\displaystyle{ n=1, }[/math] [math]\displaystyle{ 3 }[/math], 5, and [math]\displaystyle{ 7, }[/math] again with [math]\displaystyle{ \beta =2\times 10^{4}\lt math\gt and }[/math]b=50.</math> While the eigenfunctions do not have a simple shape, increasing oscillation is apparent as [math]\displaystyle{ n }[/math] increases.

Inner products

A biorthogonal system with respect to the energy inner product (\ref {energyinnerprod}) is formed by the eigenfunctions of [math]\displaystyle{ \mathcal{B} }[/math], [math]\displaystyle{ \Phi ^{+}\left( \lambda _{n},x\right) ,\lt math\gt and the eigenfunctions of }[/math]\mathcal{B} ^{\ast },[math]\displaystyle{ }[/math]\hat{\Phi}^{+}\left( \lambda _{n},x\right) </math>. To normalise the biorthogonal system, the inner product of [math]\displaystyle{ \Phi ^{+}\left( \lambda _{n},x\right) \lt math\gt and }[/math]\hat{\Phi}^{+}\left( \lambda _{n},x\right) </math> has to be determined. From the definition of the energy inner product (\ref {energyinnerprod})

[math]\displaystyle{ \begin{matrix} \left\langle \Phi \left( \lambda _{m},x\right) ,\hat{\Phi}\left( \lambda _{n},x\right) \right\rangle _{\mathcal{H}} &=&\int_{-b}^{b}\partial _{x}\phi ^{+}\left( \lambda _{m},x\right) \left( \partial _{x}\hat{\phi}^{+}\left( \lambda _{n}^{\ast },x\right) \right) ^{\ast }dx (20) \\ &&+\int_{-b}^{b}(1+P)i\zeta ^{+}\left( \lambda _{m}^{\ast },x\right) \left( i \hat{\zeta}^{+}\left( \lambda _{n}^{\ast },x\right) \right) ^{\ast }dx. \nonumber \end{matrix} }[/math]

The two integrals in (20) are considered separately. The first is

[math]\displaystyle{ \begin{matrix} &&\int_{-b}^{b}\partial _{x}\phi ^{+}\left( \lambda _{m},x\right) \left( \partial _{x}\hat{\phi}^{+}\left( \lambda _{n}^{\ast },x\right) \right) ^{\ast }dx \\ &=&\int_{-b}^{b}\left( \sum_{j=1}^{6}\mu _{j}\left( \lambda _{m}\right) \alpha _{j}e^{\mu _{j}\left( \lambda _{m}\right) x}\right) \left( \sum_{k=1}^{6}\mu _{k}\left( \lambda _{m}\right) \alpha _{k}e^{\mu _{k}\left( \lambda _{n}\right) x}\right) dx \\ &=&\sum_{j=1}^{6}\sum_{k=1}^{6}\int_{-b}^{b}-\mu _{j}\left( \lambda _{m}\right) ^{2}\alpha _{j}e^{\mu _{j}\left( \lambda _{m}\right) x}\alpha _{k}e^{\mu _{k}\left( \lambda _{n}\right) x}dx \\ &=&\sum_{j=1}^{6}\sum_{k=1}^{6}-\mu _{j}\left( \lambda _{m}\right) ^{2}\alpha _{j}\alpha _{k}\left( \frac{e^{\left( \mu _{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) \right) b}-e^{-\left( \mu _{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) \right) b} }{\mu _{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) } \right) \nonumber \end{matrix} }[/math]

and the second is

[math]\displaystyle{ \begin{matrix} &&\int_{-b}^{b}(1+P)i\zeta ^{+}\left( \lambda _{m}^{\ast },x\right) \left( i \hat{\zeta}^{+}\left( \lambda _{n}^{\ast },x\right) \right) ^{\ast }dx \\ &=&\int_{-b}^{b}\left( \sum_{j=1}^{6}-\frac{\alpha _{j}}{\lambda _{m}} \left( \mu _{j}\left( \lambda _{m}\right) ^{2}+\beta \mu _{j}\left( \lambda _{m}\right) ^{6}\right) e^{\mu _{j}\left( \lambda _{m}\right) x}\right) \left( \sum_{k=1}^{6}-\frac{\alpha _{k}}{\lambda _{n}}\mu _{k}\left( \lambda _{n}\right) ^{2}e^{\mu _{k}\left( \lambda _{n}\right) x}\right) dx \\ &=&\sum_{j=1}^{6}\sum_{k=1}^{6}\int_{-b}^{b}\frac{\alpha _{j}\alpha _{k}}{ \lambda _{m}\lambda _{n}}\left( \mu _{j}\left( \lambda _{m}\right) ^{2}+\beta \mu _{j}\left( \lambda _{m}\right) ^{6}\right) \mu _{k}\left( \lambda _{n}\right) ^{2}e^{\mu _{j}\left( \lambda _{m}\right) x}e^{\mu _{k}\left( \lambda _{n}\right) x}dx \\ &=&\sum_{j=1}^{6}\sum_{k=1}^{6}\frac{\alpha _{j}\alpha _{k}}{\lambda _{m}\lambda _{n}}\left( \mu _{j}\left( \lambda _{m}\right) ^{2}+\beta \mu _{j}\left( \lambda _{m}\right) ^{6}\right) \mu _{k}\left( \lambda _{n}\right) ^{2}\left( \frac{e^{\left( \mu _{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) \right) b}-e^{-\left( \mu _{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) \right) b}}{\mu _{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) }\right) \\ &=&\sum_{j=1}^{6}\sum_{k=1}^{6}\frac{\alpha _{j}\alpha _{k}}{\lambda _{m}\lambda _{n}}\left( -\lambda _{m}^{2}\right) \mu _{k}\left( \lambda _{n}\right) ^{2}\left( \frac{e^{\left( \mu _{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) \right) b}-e^{-\left( \mu _{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) \right) b}}{\mu _{j}\left( \lambda _{m}\right) +\mu _{k}\left( \lambda _{n}\right) }\right) . \nonumber \end{matrix} }[/math]

Therefore the calculation of the inner product in equation (\ref {innerprodall}) does not require numerical integration.

Solution

By solving (18) using the eigenfunctions of [math]\displaystyle{ \mathcal{B} }[/math] and [math]\displaystyle{ \mathcal{B}^{\ast } }[/math] we find the evolution of the plate, from an initial displacement [math]\displaystyle{ U_{0}(x), }[/math] is

[math]\displaystyle{ U\left( x,t\right) =\sum_{n=-\infty }^{\infty }e^{i\lambda _{n}t}\frac{ \left\langle U_{0}\left( x\right) ,\hat{\Phi}\left( \lambda _{n},x\right) \right\rangle _{\mathcal{H}}}{\left\langle \Phi \left( \lambda _{n},x\right) ,\hat{\Phi}\left( \lambda _{n},x\right) \right\rangle _{\mathcal{H}}}\Phi \left( \lambda _{n},x\right) . (21) }[/math]

The inner product of [math]\displaystyle{ U_{0} }[/math] with the eigenfunction [math]\displaystyle{ \hat{\Phi}\left( \lambda _{n},x\right) }[/math] is the only quantity left to compute in (\ref {evolutionB}). This inner product is written

[math]\displaystyle{ \left\langle U_{0}\left( x\right) ,\hat{\Phi}\left( \lambda _{n},x\right) \right\rangle _{\mathcal{H}}=\sum\limits_{j=1}^{6}\left( \begin{matrix}{c} \alpha _{j}\dint_{-b}^{b}-\mu _{j}\left( \lambda _{n}\right) ^{2}e^{\mu _{j}\left( \lambda _{n}\right) x}\phi _{0}\left( x\right) dx+\left. \alpha _{j}\mu _{j}\left( \lambda _{n}\right) e^{\mu _{j}\left( \lambda _{n}\right) x}\phi _{0}\left( x\right) \right| _{-b}^{b} \\ \alpha _{j}\lambda _{n}\dint_{-b}^{b}e^{\mu _{j}\left( \lambda _{n}\right) x}i\zeta _{0}\left( x\right) dx \end{matrix} \right) (22) }[/math]

and the integrals must be evaluated by numerical integration. The solutions calculated using the Lax-Phillips scattering theory are identical to those found using the self-adjoint operator method and for this reason no further figures are shown.\pagebreak

Conclusions

The spectral theory of a linear thin plate floating on shallow water has been derived. Two spectral-theory solutions have been presented which determine the time-dependent motion of the thin plate. The first method was based on self-adjoint operator theory and the second method was based on the Lax-Phillips scattering. The self-adjoint method solved both the wave forcing and the free plate problem while the Lax-Phillips method only solved for a free plate. The eigenfunctions for the self-adjoint method are orthogonal and the eigenvalues are continuous and consist of all [math]\displaystyle{ \mathbb{R} , }[/math] which makes the calculation of the eigenvalues trivial. The Lax-Phillips method has discrete eigenvalues which must be calculated numerically and the system of eigenfunctions is biorthogonal. The advantage of the Lax-Phillips method is that the modes of vibration of the plate-water system and their frequency and rate of decay are found. While the relative speeds of the two methods depends of the exact way in which they are implemented, the Lax-Phillips method should be considerably faster if the eigenvalues have been determined.

The development of a spectral theory for more complicated hydroelastic problems remains a major challenge. While this theory must be more complicated than that presented here, many features can be expected to remain. For example, ohmatsuVLFS has shown that the single frequency solutions can be used to solve certain time-dependent problems and \cite {Hazard} have shown that modes, in which the solution can be expanded, exist for other hydroelastic systems.

{\Large Acknowledgments}

\begin{acknowledgment} I would like to thank the anonymous reviewers, Dr. Kathy Ruggerio, and Prof. James Sneyd for their very helpful comments. Also, Prof. Boris Pavlov for explaining the Lax-Phillips scattering. \pagebreak \end{acknowledgment}

\bibliographystyle{jfm} \bibliography{mike,others} \pagebreak

\begin{center} {\huge Figure Captions} \end{center}

\textsc{Figure} 1. Schematic diagram of a thin plate floating on shallow water and the coordinates and dimensions of the problem.

\textsc{Figure} 2. The evolution of the displacement due to a pulse travelling to the right for the times shown. The plate occupies the region [math]\displaystyle{ -50\leq x\leq 50\lt math\gt and is shown by the bold line. }[/math]\beta =2\times 10^{4},b=50. </math>

\textsc{Figure} 3. The evolution of the potential due to a pulse travelling to the right for the times shown. The plate occupies the region [math]\displaystyle{ -50\leq x\leq 50\lt math\gt and is shown by the bold line. }[/math]\beta =2\times 10^{4},b=50.</math>

\textsc{Figure} 4. The evolution of the displacement for a plate released at [math]\displaystyle{ t=0 }[/math] for the times shown. The plate occupies the region [math]\displaystyle{ -50\leq x\leq 50 }[/math] and is shown by the bold line. [math]\displaystyle{ \beta =2\times 10^{4},b=50. }[/math]

\textsc{Figure} 5. The evolution of the displacement for a plate released at [math]\displaystyle{ t=0 }[/math] for the times shown. The plate occupies the region [math]\displaystyle{ -50\leq x\leq 50 }[/math] and is shown by the bold line. [math]\displaystyle{ \beta =2\times 10^{4},b=50. }[/math]

\textsc{Figure} 6. The location of the first 19 eigenvalues [math]\displaystyle{ \lambda _{n} }[/math] of [math]\displaystyle{ \mathcal{B} }[/math] for [math]\displaystyle{ \beta =2\times 10^{4},\;b=50. }[/math]

\textsc{Figure} 7. The real (solid) and imaginary (dashed) parts of the resonance eigenfunctions for [math]\displaystyle{ n=1, }[/math] [math]\displaystyle{ 3 }[/math], 5, and [math]\displaystyle{ 7 }[/math] as shown. [math]\displaystyle{ \beta =2\times 10^{4},\;b=50. }[/math]


\pagebreak

\FRAME{ftFU}{5.5434in}{1.6544in}{0pt}{\Qcb{}}{\Qlb{shallow}}{shallow.eps}{ \special{language "Scientific Word";type "GRAPHIC";display "ICON";valid_file "F";width 5.5434in;height 1.6544in;depth 0pt;original-width 6.9297in;original-height 1.6544in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename 'shallow.eps';file-properties "XNPEU";}}

Some nonsense\pagebreak

\FRAME{ftF}{4.6164in}{3.6296in}{0pt}{}{\Qlb{incomingspecplot1}}{ incomingspecplot1.eps}{\special{language "Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width 4.6164in;height 3.6296in;depth 0pt;original-width 6.8632in;original-height 5.3869in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename 'incomingspecplot1.eps';file-properties "XNPEU";}}

\FRAME{ftF}{4.5602in}{3.6288in}{0pt}{}{\Qlb{incomingspecplot2}}{ incomingspecplot2.eps}{\special{language "Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width 4.5602in;height 3.6288in;depth 0pt;original-width 6.7793in;original-height 5.386in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename 'incomingspecplot2.eps';file-properties "XNPEU";}}

\FRAME{fthFU}{4.6155in}{3.6357in}{0pt}{\Qcb{{}}}{\Qlb{spectral1}}{ spectral1.eps}{\special{language "Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width 4.6155in;height 3.6357in;depth 0pt;original-width 6.8493in;original-height 5.3852in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename 'spectral1.eps';file-properties "XNPEU";}}

\FRAME{fthFU}{4.5887in}{3.608in}{0pt}{\Qcb{{}}}{\Qlb{spectral2}}{ spectral2.eps}{\special{language "Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width 4.5887in;height 3.608in;depth 0pt;original-width 6.8493in;original-height 5.3852in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename 'spectral2.eps';file-properties "XNPEU";}}

\FRAME{ftFU}{4.5982in}{3.7853in}{0pt}{\Qcb{{}}}{\Qlb{spectrum}}{spectrum.eps }{\special{language "Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width 4.5982in;height 3.7853in;depth 0pt;original-width 6.8165in;original-height 5.5486in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename 'spectrum.eps';file-properties "XNPEU";}}

\FRAME{ftbpFU}{4.6164in}{3.7377in}{0pt}{\Qcb{{}}}{\Qlb{eigfunctions}}{ eigfunctions.eps}{\special{language "Scientific Word";type "GRAPHIC";maintain-aspect-ratio TRUE;display "ICON";valid_file "F";width 4.6164in;height 3.7377in;depth 0pt;original-width 6.845in;original-height 5.5365in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename 'eigfunctions.eps';file-properties "XNPEU";}}