Difference between revisions of "Kagemoto and Yue Interaction Theory"

From WikiWaves
Jump to navigationJump to search
Line 8: Line 8:
 
   
 
   
 
[[Category:Linear Water-Wave Theory]]
 
[[Category:Linear Water-Wave Theory]]
 
 
 
 
 
 
We extend the finite depth interaction theory of [[kagemoto86]] to
 
water of infinite depth and bodies of arbitrary geometry. The sum
 
over the discrete roots of the dispersion equation in the finite depth
 
theory becomes
 
an integral in the infinite depth theory. This means that the infinite
 
dimensional diffraction
 
transfer matrix
 
in the finite depth theory must be replaced by an integral
 
operator. In the numerical solution of the equations, this
 
integral operator is approximated by a sum and a linear system
 
of equations is obtained. We also show how the calculations
 
of the diffraction transfer matrix for bodies of arbitrary
 
geometry developed by [[goo90]] can be extended to
 
infinite depth, and how the diffraction transfer matrix for rotated bodies can
 
be easily calculated. This interaction theory is applied to the wave forcing
 
of multiple ice floes and a method to solve
 
the full diffraction problem in this case is presented. Convergence
 
studies comparing the interaction method with the full diffraction
 
calculations and the finite and infinite depth interaction methods are
 
carried out. 
 
 
 
 
==Introduction==
 
The scattering of water waves by floating or submerged
 
bodies is of wide practical importance.
 
Although the problem is non-linear, if the
 
wave amplitude is sufficiently small, the
 
problem can be linearised. The linear problem is still the basis of
 
the engineering design of most off-shore structures and is the standard
 
model of geophysical phenomena such as the wave forcing of ice floes.
 
While analytic solutions have been found for simplified problems
 
(especially for simple geometries or in two dimensions) the full
 
three-dimensional linear diffraction problem can only be solved by
 
numerical methods involving the discretisation of the body's surface.
 
The resulting linear system of equations has a dimension equal to the
 
number of unknowns used in the discretisation of the body.
 
 
If more than one body is present, all bodies scatter the
 
incoming waves. Therefore, the scattered wave from one body is
 
incident upon all the others and, given that they are not too far apart,
 
this notably changes the total incident wave upon them.
 
Therefore, the diffraction calculation must be conducted
 
for all bodies simultaneously. Since each body must be discretised this
 
can lead to a very large number of unknowns. However, the scattered
 
wavefield can be represented in an eigenfunction basis with a
 
comparatively small number of unknowns. If we can express the problem in this
 
basis, using what is known as an {\em Interaction Theory},
 
the number of unknowns can be much reduced, especially if there is a large
 
number of bodies.
 
 
The first interaction theory that was not based on an approximation was
 
the interaction theory
 
developed by [[kagemoto86]]. Kagemoto and Yue found an exact
 
algebraic method to solve the linear wave problem for vertically
 
non-overlapping bodies in water of finite depth.  The only restriction
 
of their theory was that the smallest escribed circle for each body must not
 
overlap any other body. The interaction of the bodies was
 
accounted for by taking the scattered wave of each body to be the
 
incident wave upon all other bodies (in addition to the ambient
 
incident wave). Furthermore, since the cylindrical eigenfunction expansions
 
are local, these were mapped from one body to another using
 
Graf's addition theorem for Bessel functions.
 
Doing this for all bodies,
 
\citeauthor{kagemoto86} were able to solve for the
 
coefficients of the scattered wavefields of all bodies simultaneously.
 
The only difficulty with this method was that the
 
solutions of the single diffraction problems had to be available in
 
the cylindrical eigenfunction expansion of an outgoing
 
wave. \citeauthor{kagemoto86} therefore only solved for axisymmetric
 
bodies because the single diffraction solution for axisymmetric
 
bodies was
 
available in the required representation. 
 
 
The extension of the Kagemoto and Yue scattering theory to bodies of
 
arbitrary geometry was performed by [[goo90]] who found a general
 
way to solve the single diffraction problem in the required
 
cylindrical eigenfunction representation. They used a
 
representation of the finite depth free surface Green's function in
 
the eigenfunction expansion of cylindrical outgoing waves
 
centred at an arbitrary point of the water surface (above the
 
body's mean centre position in this case). This Green's function was
 
presented by
 
[[black75]] and further investigated by [[fenton78]]
 
who corrected some statements about the Green's function which Black had
 
made. This Green's function is
 
based on the cylindrical
 
eigenfunction expansion of the finite depth free surface Green's
 
function given by [[john2]]. The results of
 
\citeauthor{goo90} were recently used by [[chakrabarti00]] to
 
solve for arrays of cylinders which can be divided into modules.
 
 
The development of the Kagemoto and Yue interaction theory was
 
motivated by problems
 
in off-shore engineering. However, the theory can also be applied
 
to the geophysical problem of wave scattering by ice floes. At the
 
interface of the open and frozen ocean an interfacial region known
 
as the Marginal Ice Zone (MIZ) forms. The MIZ largely controls
 
the interaction of the open and frozen ocean, especially the interaction
 
through wave processors. This is because
 
the MIZ consists of vast fields
 
of ice floes whose size is comparable to the dominant wavelength, which
 
means that it strongly scatters incoming waves. A method of
 
solving for the wave response of a single ice floe of arbitrary
 
geometry in water of infinite depth was presented by
 
[[JGR02]]. The ice floe was modelled as a floating, flexible
 
thin plate and its motion was expanded in the free plate modes of
 
vibration. Converting the problem for the water into an integral
 
equation and substituting the free modes, a system of equations for
 
the coefficients in the modal expansion was obtained. However,
 
to understand wave propagation and scattering in the MIZ we need to
 
understand the way in which large numbers of interacting
 
ice floes scatter waves. For this reason, we require an interaction
 
theory. While the Kagemoto and Yue interaction theory could be used,
 
their theory requires that the water depth is finite.
 
While the water depth in the Marginal Ice Zone varies,
 
it is generally located far from shore above the deep ocean. This means
 
that the finite depth must be chosen large in order to be able to
 
apply their theory.  Furthermore, when ocean waves propagate beneath
 
an ice floe the wavelength is increased so that it becomes more
 
difficult to make the
 
water depth sufficiently deep that it may be approximated as infinite. For this
 
reason, in this paper we develop  the equivalent interaction theory to
 
Kagemoto and Yue's in infinite depth. Also, because of
 
the complicated geometry of an ice floe, this interaction theory is
 
for bodies of arbitrary geometry. 
 
 
In the first part of this paper Kagemoto and Yue's interaction theory
 
is extended to water of infinite depth. We represent the incident and
 
scattered potentials in the cylindrical eigenfunction expansions
 
and we use an analogous infinite depth Green's function
 
to the one used by \citeauthor{goo90} \cite[given by][]{malte03}. We
 
show how the infinite
 
depth diffraction transfer matrices can be obtained with the use of
 
this Green's function and we illustrate how the rotation of a body
 
about its mean centre position in the plane can be accounted for without
 
recalculating the diffraction transfer matrix.
 
 
In the second part of the paper, using \citeauthor{JGR02}'s single
 
floe result,
 
the full diffraction calculation for the motion and scattering from many
 
interacting ice floes is calculated and presented. For two square
 
interacting ice floes the convergence of the method obtained from the
 
developed interaction theory is compared to the result of the full
 
diffraction calculation. The solutions of more than two interacting
 
ice floes and of other shapes in different arrangements are presented as well.
 
We also compare the convergence of the finite depth and infinite
 
depth methods in
 
deep water.
 
 
\section{The extension of Kagemoto and Yue's interaction
 
theory to bodies of arbitrary shape in water of infinite depth}
 
 
[[kagemoto86]] developed an interaction theory for
 
vertically non-overlapping axisymmetric structures in water of finite
 
depth. While their theory was valid for bodies of
 
arbitrary geometry, they did not develop all the necessary
 
details to apply the theory to arbitrary bodies.
 
The only requirements to apply this scattering theory is
 
that the bodies are vertically non-overlapping and
 
that the smallest cylinder which completely contains each body does not
 
intersect with any other body.
 
In this section we will extend their theory to bodies of
 
arbitrary geometry in water of infinite depth. The extension of
 
\citeauthor{kagemoto86}'s finite depth interaction theory to bodies of
 
arbitrary geometry was accomplished by [[goo90]].
 
 
 
The interaction theory begins by representing the scattered potential
 
of each body in the cylindrical eigenfunction expansion. Furthermore,
 
the incoming potential is also represented in the cylindrical
 
eigenfunction expansion. The operator which maps the incoming and
 
outgoing representation is called the diffraction transfer matrix and
 
is different for each body.
 
Since these representations are local to each body, a mapping of
 
the eigenfunction representations between different bodies
 
is required. This operator is called the coordinate transformation
 
matrix.
 
 
The cylindrical eigenfunction expansions will be introduced before we
 
derive a system of
 
equations for the coefficients of the scattered wavefields. Analogously to
 
[[kagemoto86]], we represent the scattered wavefield of
 
each body as an incoming wave upon all other bodies. The addition of
 
the ambient incident wave yields the complete incident potential and
 
with the use of diffraction transfer matrices which relate the
 
coefficients of the incident potential to those of the scattered
 
wavefield a system of equations for the unknown coefficients of the
 
scattered wavefields of all bodies is derived.
 
 
 
===Eigenfunction expansion of the potential===
 
The equations of motion for the water are derived from the linearised
 
inviscid theory. Under the assumption of irrotational motion the
 
velocity vector field of the water can be written as the gradient
 
field of a scalar velocity potential <math>\Phi</math>. Assuming that the motion
 
is time-harmonic with the radian frequency <math>\omega</math> the
 
velocity potential can be expressed as the real part of a complex
 
quantity,
 
<center><math> (time)
 
\Phi(\mathbf{y},t) = \Re \{\phi (\mathbf{y}) \mathrm{e}^{-\mathrm{i}\omega t} \}.
 
 
To simplify notation, <math>\mathbf{y} = (x,y,z)</math> will always denote a point
 
in the water, which is assumed infinitely deep, while <math>\mathbf{x}</math> will
 
always denote a point of the undisturbed water surface assumed at <math>z=0</math>.
 
 
The problem consists of <math>N</math> vertically non-overlapping bodies, denoted
 
by <math>\Delta_j</math>, which are sufficiently far apart that there is no
 
intersection of the smallest cylinder which contains each body with
 
any other body. Each body is subject to an incident wavefield which is
 
incoming, responds to this wavefield and produces a scattered wave field which
 
is outgoing. Both the incident and scattered potential corresponding
 
to these wavefields can be represented in the cylindrical
 
eigenfunction expansion valid outside of the escribed cylinder of the
 
body. Let <math>(r_j,\theta_j,z)</math> be the local cylindrical coordinates of
 
the <math>j</math>th body, <math>\Delta_j</math>, <math>j \in \{1, \ldots , N\}</math>, and <math>\alpha =
 
\omega^2/g<math> where </math>g</math> is the acceleration due to gravity. Figure
 
(fig:floe_tri) shows these coordinate systems for two bodies.
 
 
\begin{figure}
 
\begin{center}
 
\includegraphics[height=5.5cm]{floe_tri}
 
\caption{Plan view of the relation between two bodies.} (fig:floe_tri)
 
\end{center}
 
\end{figure}
 
 
The scattered potential of body <math>\Delta_j</math> can be expanded in
 
cylindrical eigenfunctions,
 
<center><math> (basisrep_out)
 
\phi_j^\mathrm{S} (r_j,\theta_j,z) &= \mathrm{e}^{\alpha z} \sum_{\nu = -
 
\infty}^{\infty} A_{0 \nu}^j H_\nu^{(1)} (\alpha r_j) \mathrm{e}^{\mathrm{i}\nu
 
\theta_j}\\
 
&\quad + \int\limits_0^{\infty} \big( \cos \eta z + \frac{\alpha}{\eta}
 
\sin \eta z \big) \sum_{\nu = -
 
\infty}^{\infty} A_{\nu}^j (\eta) K_\nu (\eta r_j) \mathrm{e}^{\mathrm{i}\nu
 
\theta_j} \mathrm{d}\eta,
 
 
where the coefficients <math>A_{0 \nu}^j</math> for the propagating modes are
 
discrete and the coefficients <math>A_{\nu}^j (\cdot)</math> for the decaying
 
modes are functions. <math>H_\nu^{(1)}</math> and <math>K_\nu</math> are the Hankel function
 
of the first kind and the modified Bessel function of the second kind
 
respectively, both of order <math>\nu</math> as defined in [[abr_ste]].
 
The incident potential upon body <math>\Delta_j</math> can be expanded in
 
cylindrical eigenfunctions,
 
<center><math> (basisrep_in)
 
\phi_j^\mathrm{I} (r_j,\theta_j,z) &= \mathrm{e}^{\alpha z} \sum_{\mu = -
 
\infty}^{\infty} D_{0 \mu}^j J_\mu (\alpha r_j) \mathrm{e}^{\mathrm{i}\mu
 
\theta_j}\\
 
& \quad + \int\limits_0^{\infty} \big( \cos \eta z + \frac{\alpha}{\eta}
 
\sin \eta z \big) \sum_{\mu = -
 
\infty}^{\infty} D_{\mu}^j (\eta) I_\mu (\eta r_j) \mathrm{e}^{\mathrm{i}\mu
 
\theta_j} \mathrm{d}\eta,
 
 
where the coefficients <math>D_{0 \mu}^j</math> for the propagating modes are
 
discrete and the coefficients <math>D_{\mu}^j (\cdot)</math> for the decaying
 
modes are functions. <math>J_\mu</math> and <math>I_\mu</math> are the Bessel function and
 
the modified Bessel function respectively, both of the first kind and
 
order <math>\mu</math>. To simplify the notation, from now on <math>\psi(z,\eta)</math> will
 
denote the vertical eigenfunctions corresponding to the decaying modes,
 
<center><math>
 
\psi(z,\eta) = \cos \eta z + \alpha / \eta \, \sin \eta z.
 
 
 
===The interaction in water of infinite depth===
 
Following the ideas of [[kagemoto86]], a system of equations for the unknown
 
coefficients and coefficient functions of the scattered wavefields
 
will be developed. This system of equations is based on transforming the
 
scattered potential of <math>\Delta_j</math> into an incident potential upon
 
<math>\Delta_l</math> (<math>j \neq l</math>). Doing this for all bodies simultaneously,
 
and relating the incident and scattered potential for each body, a system
 
of equations for the unknown coefficients will be developed.
 
 
The scattered potential <math>\phi_j^{\mathrm{S}}</math> of body <math>\Delta_j</math> needs to be
 
represented in terms of the incident potential <math>\phi_l^{\mathrm{I}}</math>
 
upon <math>\Delta_l</math>, <math>j \neq l</math>. From figure
 
(fig:floe_tri) we can see that this can be accomplished by using
 
Graf's addition theorem for Bessel functions given in
 
\citet[eq. 9.1.79]{abr_ste},
 
\begin{subequations} (transf)
 
<center><math>\begin{matrix} (transf_h)
 
H_\nu^{(1)}(\alpha r_j) \mathrm{e}^{\mathrm{i}\nu (\theta_j - \vartheta_{jl})} &=
 
\sum_{\mu = - \infty}^{\infty} H^{(1)}_{\nu + \mu} (\alpha R_{jl}) \,
 
J_\mu (\alpha r_l) \mathrm{e}^{\mathrm{i}\mu (\pi - \theta_l + \vartheta_{jl})},
 
\quad j \neq l,\\
 
(transf_k)
 
K_\nu(\eta r_j) \mathrm{e}^{\mathrm{i}\nu (\theta_j - \vartheta_{jl})} &= \sum_{\mu = -
 
\infty}^{\infty} K_{\nu + \mu} (\eta R_{jl}) \, I_\mu (\eta r_l) \mathrm{e}^{\mathrm{i}\mu
 
(\pi - \theta_l + \vartheta_{jl})}, \quad j \neq l,
 
\end{matrix}</math></center>
 
\end{subequations}
 
which is valid provided that <math>r_l < R_{jl}</math>. This limitation
 
only requires that the escribed cylinder of each body <math>\Delta_l</math> does
 
not enclose any other origin <math>O_j</math> (<math>j \neq l</math>). However, the
 
expansion of the scattered and incident potential in cylindrical
 
eigenfunctions is only valid outside the escribed cylinder of each
 
body. Therefore the condition that the
 
escribed cylinder of each body <math>\Delta_l</math> does not enclose any other
 
origin <math>O_j</math> (<math>j \neq l</math>) is superseded by the more rigorous
 
restriction that the escribed cylinder of each body may not contain any
 
other body. Making use of the equations  (transf)
 
the scattered potential of <math>\Delta_j</math> can be expressed in terms of the
 
incident potential upon <math>\Delta_l</math>,
 
<center><math>\begin{matrix}
 
\phi_j^{\mathrm{S}} (r_l,\theta_l,z) &=
 
\mathrm{e}^{\alpha z} \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j
 
\sum_{\mu = -\infty}^{\infty} H_{\nu - \mu}^{(1)} (\alpha R_{jl})
 
J_\mu (\alpha r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l} \mathrm{e}^{\mathrm{i}(\nu-\mu)
 
\vartheta_{jl}}\\
 
& \quad + \int\limits_0^{\infty} \psi (z,\eta) \sum_{\nu = -
 
\infty}^{\infty} A_{\nu}^j (\eta) \sum_{\mu = -\infty}^{\infty}
 
(-1)^\mu K_{\nu-\mu} (\eta R_{jl}) I_\mu (\eta r_l) \mathrm{e}^{\mathrm{i}\mu
 
\theta_l}  \mathrm{e}^{\mathrm{i}(\nu-\mu) \vartheta_{jl}} \mathrm{d}\eta\\
 
&= \mathrm{e}^{\alpha z} \sum_{\mu = -\infty}^{\infty} \Big[ \sum_{\nu = -
 
\infty}^{\infty} A_{0\nu}^j H_{\nu - \mu}^{(1)} (\alpha R_{jl}) \mathrm{e}^{\i
 
(\nu-\mu) \vartheta_{jl}} \Big] J_\mu (\alpha r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l}\\
 
& + \int\limits_0^{\infty} \psi (z,\eta) \sum_{\mu = -\infty}^{\infty}
 
\Big[ \sum_{\nu = - \infty}^{\infty} A_{\nu}^j (\eta) (-1)^\mu K_{\nu-\mu}
 
(\eta R_{jl}) \mathrm{e}^{\mathrm{i}(\nu-\mu)
 
\vartheta_{jl}} \Big] I_\mu (\eta r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l} \mathrm{d}\eta.
 
\end{matrix}</math></center>
 
The ambient incident wavefield <math>\phi^{\mathrm{In}}</math> can also be
 
expanded in the eigenfunctions corresponding to the incident wavefield upon
 
<math>\Delta_l</math>. A detailed illustration of how to accomplish this will be given
 
later. Let <math>D_{l0\mu}^{\mathrm{In}}</math> denote the coefficients of this
 
ambient incident wavefield corresponding to the propagating modes and
 
<math>D_{l\mu}^{\mathrm{In}} (\cdot)</math>  denote the coefficients functions
 
corresponding to the decaying modes (which are identically zero) of
 
the incoming eigenfunction expansion for <math>\Delta_l</math>. The total
 
incident wavefield upon body <math>\Delta_j</math> can now be expressed as
 
<center><math>\begin{matrix}
 
&\phi_l^{\mathrm{I}}(r_l,\theta_l,z) = \phi^{\mathrm{In}}(r_l,\theta_l,z) +
 
\sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \, \phi_j^{\mathrm{S}}
 
(r_l,\theta_l,z)\\
 
&= \mathrm{e}^{\alpha z} \sum_{\mu = -\infty}^{\infty} \Big[
 
D_{l0\mu}^{\mathrm{In}} + \sum_{\genfrac{}{}{0pt}{}{j=1}{j
 
\neq l}}^{N} \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j
 
H_{\nu-\mu}^{(1)} (\alpha R_{jl}) \mathrm{e}^{\i
 
(\nu - \mu) \vartheta_{jl}} \Big] J_\mu (\alpha r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l}\\
 
& + \int\limits_0^{\infty} \psi (z,\eta) \sum_{\mu =
 
-\infty}^{\infty} \Big[  D_{l\mu}^{\mathrm{In}}(\eta) +
 
\sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq  l}}^{N} \sum_{\nu =
 
-\infty}^{\infty} A_{\nu}^j (\eta) (-1)^\mu K_{\nu - \mu}  (\eta
 
R_{jl}) \mathrm{e}^{\mathrm{i}(\nu - \mu) \vartheta_{jl}} \Big] I_\mu (\eta r_l)
 
\mathrm{e}^{\mathrm{i}\mu \theta_l} \mathrm{d}\eta.
 
\end{matrix}</math></center>
 
The coefficients of the total incident potential upon <math>\Delta_l</math> are
 
therefore given by
 
\begin{subequations} (inc_coeff)
 
<center><math>\begin{matrix}
 
D_{0\mu}^l &= D_{l0\mu}^{\mathrm{In}} + \sum_{\genfrac{}{}{0pt}{}{j=1}{j
 
\neq l}}^{N} \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j
 
H_{\nu-\mu}^{(1)} (\alpha R_{jl}) \mathrm{e}^{\i
 
(\nu - \mu) \vartheta_{jl}},\\
 
D_{\mu}^l(\eta) &= D_{l\mu}^{\mathrm{In}}(\eta) +
 
\sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq  l}}^{N} \sum_{\nu =
 
-\infty}^{\infty} A_{\nu}^j (\eta) (-1)^\mu K_{\nu - \mu}  (\eta
 
R_{jl}) \mathrm{e}^{\mathrm{i}(\nu - \mu) \vartheta_{jl}}.
 
\end{matrix}</math></center>
 
\end{subequations}
 
 
In general, it is possible to relate the total incident and scattered
 
partial waves for any body through the diffraction characteristics of
 
that body in isolation. There exist diffraction transfer operators
 
<math>B_l</math> that relate the coefficients of the incident and scattered
 
partial waves, such that
 
<center><math> (eq_B)
 
A_l = B_l (D_l), \quad l=1, \ldots, N,
 
 
where <math>A_l</math> are the scattered modes due to the incident modes <math>D_l</math>.
 
In the case of a countable number of modes, (i.e. when
 
the depth is finite), <math>B_l</math> is an infinite dimensional matrix. When
 
the modes are functions of a continuous variable (i.e. infinite
 
depth), <math>B_l</math> is the kernel of an integral operator.
 
For the propagating and the decaying modes respectively, the scattered
 
potential can be related by diffraction transfer operators acting in the
 
following ways,
 
\begin{subequations} (diff_op)
 
<center><math>\begin{matrix}
 
A_{0\nu}^l &= \sum_{\mu = -\infty}^{\infty} B_{l\nu\mu}^\mathrm{pp} D_{0\mu}^l
 
+ \int\limits_{0}^{\infty} \sum_{\mu = -\infty}^{\infty}
 
B_{l\nu\mu}^\mathrm{pd} (\xi) D_{\mu}^l (\xi) \mathrm{d}\xi,\\
 
A_\nu^l (\eta) &= \sum_{\mu = -\infty}^{\infty}
 
B_{l\nu\mu}^\mathrm{dp} (\eta) D_{0\mu}^l + \int\limits_{0}^{\infty}
 
\sum_{\mu = -\infty}^{\infty} B_{l\nu\mu}^\mathrm{dd} (\eta;\xi)
 
D_{\mu}^l (\xi) \mathrm{d}\xi.
 
\end{matrix}</math></center>
 
\end{subequations}
 
The superscripts <math>\mathrm{p}</math> and <math>\mathrm{d}</math> are used to distinguish
 
between propagating and decaying modes, the first superscript denotes the kind
 
of scattered mode, the second one the kind of incident mode.
 
If the diffraction transfer operators are known (their calculation
 
will be discussed later), the substitution of
 
equations  (inc_coeff) into equations  (diff_op) give the
 
required equations to determine the coefficients and coefficient
 
functions of the scattered wavefields of all bodies,
 
\begin{subequations} (eq_op)
 
<center><math>\begin{matrix}
 
&\begin{aligned}
 
&A_{0n}^l = \sum_{\mu = -\infty}^{\infty} B_{ln\mu}^\mathrm{pp}
 
\Big[ D_{l0\mu}^{\mathrm{In}} + \sum_{\genfrac{}{}{0pt}{}{j=1}{j
 
\neq l}}^{N} \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j
 
H_{\nu-\mu}^{(1)} (\alpha R_{jl}) \mathrm{e}^{\i
 
(\nu - \mu) \vartheta_{jl}} \Big]\\
 
& \ + \int\limits_{0}^{\infty} \sum_{\mu = -\infty}^{\infty}
 
B_{ln\mu}^\mathrm{pd} (\xi) \Big[D_{l\mu}^{\mathrm{In}}(\eta) +
 
\sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq  l}}^{N} \sum_{\nu =
 
-\infty}^{\infty} A_{\nu}^j (\eta) (-1)^\mu K_{\nu - \mu}  (\eta
 
R_{jl}) \mathrm{e}^{\mathrm{i}(\nu - \mu) \vartheta_{jl}} \Big] \mathrm{d}\xi,
 
\end{aligned}\\
 
&\begin{aligned}
 
&A_n^l (\eta) = \sum_{\mu = -\infty}^{\infty}
 
B_{ln\mu}^\mathrm{dp} (\eta) \Big[
 
D_{l0\mu}^{\mathrm{In}} + \sum_{\genfrac{}{}{0pt}{}{j=1}{j
 
\neq l}}^{N} \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j
 
H_{\nu-\mu}^{(1)} (\alpha R_{jl}) \mathrm{e}^{\i
 
(\nu - \mu) \vartheta_{jl}}\Big]\\
 
& \ + \int\limits_{0}^{\infty}
 
\sum_{\mu = -\infty}^{\infty} B_{ln\mu}^\mathrm{dd} (\eta;\xi)
 
\Big[ D_{l\mu}^{\mathrm{In}}(\eta) +
 
\sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq  l}}^{N} \sum_{\nu =
 
-\infty}^{\infty} A_{\nu}^j (\eta) (-1)^\mu K_{\nu - \mu}  (\eta
 
R_{jl}) \mathrm{e}^{\mathrm{i}(\nu - \mu) \vartheta_{jl}}\Big] \mathrm{d}\xi,
 
\end{aligned}
 
\end{matrix}</math></center>
 
\end{subequations}
 
<math>n \in \mathds{Z},\, l = 1, \ldots, N</math>. It has to be noted that all
 
equations are coupled so that it is necessary to solve for all
 
scattered coefficients and coefficient functions simultaneously.
 
 
For numerical calculations, the infinite sums have to be truncated and
 
the integrals must be discretised. Implying a suitable truncation, the
 
four different diffraction transfer operators can be represented by
 
matrices which can be assembled in a big matrix <math>\mathbf{B}_l</math>,
 
<center><math>
 
\mathbf{B}_l = \left[
 
\begin{matrix}{cc} \mathbf{B}_l^{\mathrm{pp}} & \mathbf{B}_l^{\mathrm{pd}}\\
 
\mathbf{B}_l^{\mathrm{dp}} & \mathbf{B}_l^{\mathrm{dd}}
 
\end{matrix} \right],
 
 
the infinite depth diffraction transfer matrix.
 
Truncating the coefficients accordingly, defining <math>{\bf a}^l</math> to be the
 
vector of the coefficients of the scattered potential of body
 
<math>\Delta_l</math>, <math>\mathbf{d}_l^{\mathrm{In}}</math> to be the vector of
 
coefficients of the ambient wavefield, and making use of a coordinate
 
transformation matrix <math>{\bf T}_{jl}</math> given by
 
\begin{subequations} (T_elem_deep)
 
<center><math>\begin{matrix}
 
({\bf T}_{jl})_{pq} &= H_{p-q}^{(1)}(\alpha R_{jl}) \, \mathrm{e}^{\mathrm{i}(p-q)
 
\vartheta_{jl}}\\
 
=for the propagating modes, and=
 
({\bf T}_{jl})_{pq} &= (-1)^{q} \, K_{p-q} (\eta R_{jl}) \, \mathrm{e}^{\i
 
(p-q) \vartheta_{jl}}
 
\end{matrix}</math></center>
 
\end{subequations}
 
for the decaying modes, a linear system of equations
 
for the unknown coefficients follows from equations  (eq_op),
 
<center><math> (eq_B_inf)
 
{\bf a}_l = {\bf \hat{B}}_l \Big( {\bf d}_l^{\mathrm{In}} +
 
\sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \trans {\bf T}_{jl} \,
 
{\bf a}_j \Big), \quad  l=1, \ldots, N,
 
 
where the left superscript <math>\mathrm{t}</math> indicates transposition.
 
The matrix <math>{\bf \hat{B}}_l</math> denotes the infinite depth diffraction
 
transfer matrix <math>{\bf B}_l</math> in which the elements associated with
 
decaying scattered modes have been multiplied with the appropriate
 
integration weights depending on the discretisation of the continuous variable.
 
 
 
 
\subsection{Calculation of the diffraction transfer matrix for bodies
 
of arbitrary geometry}
 
 
Before we can apply the interaction theory we require the diffraction
 
transfer matrices <math>\mathbf{B}_j</math> which relate the incident and the
 
scattered potential for a body <math>\Delta_j</math> in isolation.
 
The elements of the diffraction transfer matrix, <math>({\bf B}_j)_{pq}</math>,
 
are the coefficients of the
 
<math>p</math>th partial wave of the scattered potential due to a single
 
unit-amplitude incident wave of mode <math>q</math> upon <math>\Delta_j</math>.
 
 
While \citeauthor{kagemoto86}'s interaction theory was valid for
 
bodies of arbitrary shape, they did not explain how to actually obtain the
 
diffraction transfer matrices for bodies which did not have an axisymmetric
 
geometry. This step was performed by [[goo90]] who came up with an
 
explicit method to calculate the diffraction transfer matrices for bodies of
 
arbitrary geometry in the case of finite depth. Utilising a Green's
 
function they used the standard
 
method of transforming the single diffraction boundary-value problem
 
to an integral equation for the source strength distribution function
 
over the immersed surface of the body.
 
However, the representation of the scattered potential which
 
is obtained using this method is not automatically given in the
 
cylindrical eigenfunction
 
expansion. To obtain such cylindrical eigenfunction expansions of the
 
potential [[goo90]] used the representation of the free surface
 
finite depth Green's function given by [[black75]] and
 
[[fenton78]].  \citeauthor{black75} and
 
\citeauthor{fenton78}'s representation of the Green's function was based
 
on applying Graf's addition theorem to the eigenfunction
 
representation of the free surface finite depth Green's function given
 
by [[john2]]. Their representation allowed the scattered potential to be
 
represented in the eigenfunction expansion with the cylindrical
 
coordinate system fixed at the point of the water surface above the
 
mean centre position of the body.
 
 
It should be noted that, instead of using the source strength distribution
 
function, it is also possible to consider an integral equation for the
 
total potential and calculate the elements of the diffraction transfer
 
matrix from the solution of this integral equation.
 
An outline of this method for water of finite
 
depth is given by [[kashiwagi00]]. We will present
 
here a derivation of the diffraction transfer matrices for the case
 
infinite depth based on a solution
 
for the source strength distribution function. However,
 
an equivalent derivation would be possible based on the solution
 
for the total velocity potential.
 
 
To calculate the diffraction transfer matrix in infinite depth, we
 
require the representation of the infinite depth free surface Green's
 
function in cylindrical eigenfunctions,
 
<center><math> (green_inf)
 
G(r,\theta,z;s,\varphi,c) &= \frac{\mathrm{i}\alpha}{2} \,  \mathrm{e}^{\alpha (z+c)}
 
\sum_{\nu=-\infty}^{\infty} H_\nu^{(1)}(\alpha r) J_\nu(\alpha s) \mathrm{e}^{\mathrm{i}\nu
 
(\theta - \varphi)}\\ &\quad + \frac{1}{\pi^2} \int\limits_0^{\infty}
 
\psi(z,\eta) \frac{\eta^2}{\eta^2+\alpha^2} \psi(c,\eta)
 
\sum_{\nu=-\infty}^{\infty} K_\nu(\eta r) I_\nu(\eta s) \mathrm{e}^{\mathrm{i}\nu
 
(\theta - \varphi)} \mathrm{d}\eta,
 
 
<math>r > s</math>, given by [[malte03]].
 
 
We assume that we have represented the scattered potential in terms of
 
the source strength distribution <math>\varsigma^j</math> so that the scattered
 
potential can be written as
 
<center><math> (int_eq_1)
 
\phi_j^\mathrm{S}(\mathbf{y}) = \int\limits_{\Gamma_j} G
 
(\mathbf{y},\mathbf{\zeta}) \, \varsigma^j (\mathbf{\zeta})
 
\mathrm{d}\sigma_\mathbf{\zeta}, \quad \mathbf{y} \in D,
 
 
where <math>D</math> is the volume occupied by the water and <math>\Gamma_j</math> is the
 
immersed surface of body <math>\Delta_j</math>. The source strength distribution
 
function <math>\varsigma^j</math> can be found by solving an
 
integral equation. The integral equation is described in
 
[[Weh_Lait]] and numerical methods for its solution are outlined in
 
[[Sarp_Isa]].
 
Substituting the eigenfunction expansion of the Green's function
 
(green_inf) into  (int_eq_1), the scattered potential can
 
be written as
 
<center><math>\begin{matrix}
 
&\phi_j^\mathrm{S}(r_j,\theta_j,z) = \mathrm{e}^{\alpha z} \sum_{\nu = -
 
\infty}^{\infty} \bigg[ \frac{\mathrm{i}\alpha}{2}
 
\int\limits_{\Gamma_j} \mathrm{e}^{\alpha c} J_\nu(\alpha s) \mathrm{e}^{-\mathrm{i}\nu
 
\varphi} \varsigma^j(\mathbf{\zeta})
 
\mathrm{d}\sigma_\mathbf{\zeta} \bigg] H_\nu^{(1)} (\alpha r_j) \mathrm{e}^{\mathrm{i}\nu \theta_j}\\
 
& \, + \int\limits_{0}^{\infty} \psi(z,\eta) \sum_{\nu = -
 
\infty}^{\infty}  \bigg[ \frac{1}{\pi^2} \frac{\eta^2
 
}{\eta^2 + \alpha^2} \int\limits_{\Gamma_j} \psi(c,\eta) I_\nu(\eta s)
 
\mathrm{e}^{-\mathrm{i}\nu \varphi} \varsigma^j({\bf{\zeta}})
 
\mathrm{d}\sigma_{\bf{\zeta}} \bigg] K_\nu (\eta r_j) \mathrm{e}^{\mathrm{i}\nu \theta_j} \mathrm{d}\eta,
 
\end{matrix}</math></center>
 
where
 
<math>\mathbf{\zeta}=(s,\varphi,c)</math> and <math>r>s</math>.
 
This restriction implies that the eigenfunction expansion is only valid
 
outside the escribed cylinder of the body.
 
 
The columns of the diffraction transfer matrix are the coefficients of
 
the eigenfunction expansion of the scattered wavefield due to the
 
different incident modes of unit-amplitude. The elements of the
 
diffraction transfer matrix of a body of arbitrary shape are therefore given by
 
\begin{subequations} (B_elem)
 
<center><math>\begin{matrix}
 
({\bf B}_j)_{pq} &= \frac{\mathrm{i}\alpha}{2} \int\limits_{\Gamma_j}
 
\mathrm{e}^{\alpha c} J_p(\alpha s) \mathrm{e}^{-\mathrm{i}p \varphi} \varsigma_q^j(\mathbf{\zeta})
 
\mathrm{d}\sigma_\mathbf{\zeta}\\
 
=and= 
 
({\bf B}_j)_{pq} &= \frac{1}{\pi^2} \frac{\eta^2}{\eta^2 +
 
\alpha^2} \int\limits_{\Gamma_j} \psi(c,\eta) I_p(\eta s) \mathrm{e}^{-\mathrm{i}p
 
\varphi} \varsigma_q^j(\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta}
 
\end{matrix}</math></center>
 
\end{subequations}
 
for the propagating and the decaying modes respectively, where
 
<math>\varsigma_q^j(\mathbf{\zeta})</math> is the source strength distribution
 
due to an incident potential of mode <math>q</math> of the form
 
\begin{subequations} (test_modes_inf)
 
<center><math>\begin{matrix}
 
\phi_q^{\mathrm{I}}(s,\varphi,c) &=  \mathrm{e}^{\alpha c} H_q^{(1)} (\alpha
 
s) \mathrm{e}^{\mathrm{i}q \varphi}\\
 
=for the propagating modes, and= 
 
\phi_q^{\mathrm{I}}(s,\varphi,c) &= \psi(c,\eta) K_q (\eta s) \mathrm{e}^{\mathrm{i}q \varphi}
 
\end{matrix}</math></center>
 
\end{subequations}
 
for the decaying modes.
 
 
===The diffraction transfer matrix of rotated bodies===
 
 
For a non-axisymmetric body, a rotation about the mean
 
centre position in the <math>(x,y)</math>-plane will result in a
 
different diffraction transfer matrix. We will show how the
 
diffraction transfer matrix of a body rotated by an angle <math>\beta</math> can
 
be easily calculated from the diffraction transfer matrix of the
 
non-rotated body. The rotation of the body influences the form of the
 
elements of the diffraction transfer matrices in two ways. Firstly, the
 
angular dependence in the integral over the immersed surface of the
 
body is altered and, secondly, the source strength distribution
 
function is different if the body is rotated. However, the source
 
strength distribution function of the rotated body can be obtained by
 
calculating the response of the non-rotated body due to rotated
 
incident potentials. It will be shown that the additional angular
 
dependence can be easily factored out of the elements of the
 
diffraction transfer matrix.
 
 
The additional angular dependence caused by the rotation of the
 
incident potential can be factored out of the normal derivative of the
 
incident potential such that
 
<center><math>
 
\frac{\partial \phi_{q\beta}^{\mathrm{I}}}{\partial n} =
 
\frac{\partial \phi_{q}^{\mathrm{I}}}{\partial n}
 
\mathrm{e}^{\mathrm{i}q \beta},
 
 
where <math>\phi_{q\beta}^{\mathrm{I}}</math> is the rotated incident potential.
 
Since the integral equation for the determination of the source
 
strength distribution function is linear, the source strength
 
distribution function due to the rotated incident potential is thus just
 
given by
 
<center><math>
 
\varsigma_{q\beta}^j = \varsigma_q^j \, \mathrm{e}^{\mathrm{i}q \beta}.
 
 
This is also the source strength distribution function of the rotated
 
body due to the standard incident modes.
 
 
The elements of the diffraction transfer matrix <math>\mathbf{B}_j</math> are
 
given by equations  (B_elem). Keeping in mind that the body is
 
rotated by the angle <math>\beta</math>, the elements of the diffraction transfer
 
matrix of the rotated body are given by
 
\begin{subequations} (B_elem_rot)
 
<center><math>\begin{matrix}
 
(\mathbf{B}_j^\beta)_{pq} &= \frac{\mathrm{i}\alpha}{2} \int\limits_{\Gamma_j}
 
\mathrm{e}^{\alpha c} J_p(\alpha s) \mathrm{e}^{-\mathrm{i}p (\varphi+\beta)}
 
\varsigma_{q\beta}^j(\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta},\\
 
=and= 
 
(\mathbf{B}_j^\beta)_{pq} &= \frac{1}{\pi^2} \frac{\eta^2}{\eta^2 +
 
\alpha^2} \int\limits_{\Gamma_j} \psi(c,\eta) I_p(\eta s) \mathrm{e}^{-\mathrm{i}p
 
(\varphi+\beta)} \varsigma_{q\beta}^j(\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta},
 
\end{matrix}</math></center>
 
\end{subequations}
 
for the propagating and decaying modes respectively.
 
 
Thus the additional angular dependence caused by the rotation of
 
the body can be factored out of the elements of the diffraction
 
transfer matrix. The elements of the diffraction transfer matrix
 
corresponding to the body rotated by the angle <math>\beta</math>,
 
<math>\mathbf{B}_j^\beta</math>, are given by
 
<center><math> (B_rot)
 
(\mathbf{B}_j^\beta)_{pq} = (\mathbf{B}_j)_{pq} \, \mathrm{e}^{\mathrm{i}(q-p) \beta}.
 
 
As before, <math>(\mathbf{B})_{pq}</math> is understood to be the element of
 
<math>\mathbf{B}</math> which corresponds to the coefficient of the <math>p</math>th scattered
 
mode due to a unit-amplitude incident wave of mode <math>q</math>. Equation
 
(B_rot) applies to propagating and decaying modes likewise.
 
 
 
\subsection{Representation of the ambient wavefield in the eigenfunction
 
representation}
 
In Cartesian coordinates centred at the origin, the ambient wavefield is
 
given by
 
<center><math>
 
\phi^{\mathrm{In}}(x,y,z) = A \frac{g}{\omega} \, \mathrm{e}^{\mathrm{i}\alpha (x
 
\cos \chi + y \sin \chi)+ \alpha z},
 
 
where <math>A</math> is the amplitude (in displacement) and <math>\chi</math> is the
 
angle between the <math>x</math>-axis and the direction in which the wavefield travels.
 
The interaction theory requires that the ambient wavefield, which is
 
incident upon
 
all bodies, is represented in the eigenfunction expansion of an
 
incoming wave in the local coordinates of the body. The ambient wave
 
can be represented in an eigenfunction expansion centred at the origin
 
as
 
<center><math>
 
\phi^{\mathrm{In}}(x,y,z) = A \frac{g}{\omega} \mathrm{e}^{\alpha z}
 
\sum_{\mu = -\infty}^{\infty} \mathrm{e}^{\mathrm{i}\mu (\pi/2 - \theta + \chi)}
 
J_\mu(\alpha r)
 
 
\cite[p. 169]{linton01}.
 
Since the local coordinates of the bodies are centred at their mean
 
centre positions <math>O_l = (O_x^l,O_y^l)</math>, a phase factor has to be defined
 
which accounts for the position from the origin. Including this phase
 
factor the ambient wavefield at the <math>l</math>th body is given
 
by
 
<center><math>
 
\phi^{\mathrm{In}}(r_l,\theta_l,z) = A \frac{g}{\omega} \, e^{\mathrm{i}\alpha (O_x^l
 
\cos \chi + O_x^l \sin \chi)} \, \mathrm{e}^{\alpha z}
 
\sum_{\mu = -\infty}^{\infty} \mathrm{e}^{\mathrm{i}\mu (\pi/2 -  \chi)}
 
J_\mu(\alpha r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l}.
 
 
 
===Solving the resulting system of equations===
 
After the coefficient vector of the ambient incident wavefield, the
 
diffraction transfer matrices and the coordinate
 
transformation matrices have been calculated, the system of
 
equations  (eq_B_inf),
 
has to be solved. This system can be represented by the following
 
matrix equation,
 
<center><math>
 
\left[ \begin{matrix}{c}
 
{\bf a}_1\\ {\bf a}_2\\ \\ \vdots \\ \\ {\bf a}_N
 
\end{matrix} \right]
 
= \left[ \begin{matrix}{c}
 
\hat{{\bf B}}_1 {\bf d}_1^\mathrm{In}\\ \hat{{\bf B}}_2 {\bf
 
d}_2^\mathrm{In}\\ \\ \vdots \\ \\ \hat{{\bf B}}_N {\bf d}_N^\mathrm{In}
 
\end{matrix} \right]+
 
\left[ \begin{matrix}{ccccc}
 
\mathbf{0} & \hat{{\bf B}}_1 \trans {\bf T}_{21} & \hat{{\bf B}}_1
 
\trans {\bf T}_{31} & \dots & \hat{{\bf B}}_1 \trans {\bf T}_{N1}\\
 
\hat{{\bf B}}_2 \trans {\bf T}_{12} & \mathbf{0} & \hat{{\bf B}}_2
 
\trans {\bf T}_{32} & \dots & \hat{{\bf B}}_2 \trans {\bf T}_{N2}\\
 
& & \mathbf{0} & &\\
 
\vdots & & & \ddots & \vdots\\
 
& & & & \\
 
\hat{{\bf B}}_N \trans {\bf T}_{1N} & & \dots & 
 
& \mathbf{0}
 
\end{matrix} \right]
 
\left[ \begin{matrix}{c}
 
{\bf a}_1\\ {\bf a}_2\\ \\ \vdots \\ \\ {\bf a}_N
 
\end{matrix} \right],
 
 
where <math>\mathbf{0}</math> denotes the zero-matrix which is of the same
 
dimension as <math>\hat{{\bf B}}_j</math>, say <math>n</math>. This matrix equation can be
 
easily transformed into a classical <math>(N \, n)</math>-dimensional linear system of
 
equations.
 
 
 
==Finite Depth Interaction Theory==
 
 
We will compare the performance of the infinite depth interaction theory
 
with the equivalent theory for finite
 
depth. As we have stated previously, the finite depth theory was
 
developed by [[kagemoto86]] and extended to bodies of arbitrary
 
geometry by [[goo90]]. We will briefly present this theory in
 
our notation and the comparisons will be made in a later section.
 
 
In water of constant finite depth <math>d</math>, the scattered potential of a body
 
<math>\Delta_j</math> can be expanded in cylindrical eigenfunctions,
 
<center><math> (basisrep_out_d)
 
\phi_j^\mathrm{S} (r_j,\theta_j,z) &= \frac{\cosh k(z+d)}{\cosh kd}
 
\sum_{\nu = - \infty}^{\infty} A_{0 \nu}^j H_\nu^{(1)} (k r_j) \mathrm{e}^{\mathrm{i}\nu
 
\theta_j}\\
 
&\quad + \sum_{m=1}^{\infty} \frac{\cos k_m (z+d)}{\cos kd} \sum_{\nu = -
 
\infty}^{\infty} A_{m \nu}^j K_\nu (k_m r_j) \mathrm{e}^{\mathrm{i}\nu
 
\theta_j},
 
 
with discrete coefficients <math>A_{m \nu}^j</math>. The positive wavenumber <math>k</math>
 
is related to <math>\alpha</math> by the dispersion relation
 
<center><math> (eq_k)
 
\alpha = k \tanh k d,
 
 
and the values of <math>k_m</math>, <math>m>0</math>, are given as positive real roots of
 
the dispersion relation
 
<center><math> (eq_k_m)
 
\alpha + k_m \tan k_m d = 0.
 
 
The incident potential upon body <math>\Delta_j</math> can be also be expanded in
 
cylindrical eigenfunctions,
 
<center><math> (basisrep_in_d)
 
\phi_j^\mathrm{I} (r_j,\theta_j,z) &= \frac{\cosh k(z+d)}{\cosh kd}
 
\sum_{\mu = - \infty}^{\infty} D_{0 \mu}^j J_\mu (k r_j) \mathrm{e}^{\mathrm{i}\mu
 
\theta_j}\\
 
& \quad + \sum_{m=1}^{\infty} \frac{\cos k_m (z+d)}{\cos kd} \sum_{\mu = -
 
\infty}^{\infty} D_{m\mu}^j I_\mu (k_m r_j) \mathrm{e}^{\mathrm{i}\mu
 
\theta_j},
 
 
with discrete coefficients <math>D_{m\mu}^j</math>. A system of equations for the
 
coefficients of the scattered wavefields for the bodies are derived
 
in an analogous way to the infinite depth case. The derivation is
 
simpler because all the coefficients are discrete and the
 
diffraction transfer operator can be represented by an
 
infinite dimensional matrix.
 
Truncating the infinite dimensional matrix as well as the
 
coefficient vectors appropriately, the resulting system of
 
equations is given by 
 
<center><math>
 
{\bf a}_l = {\bf B}_l \Big( {\bf d}_l^\mathrm{In} +
 
\sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \trans {\bf T}_{jl} \,
 
{\bf a}_j \Big), \quad  l=1, \ldots, N,
 
 
where <math>{\bf a}_l</math> is the coefficient vector of the scattered
 
wave, <math>{\bf d}_l^\mathrm{In}</math> is the coefficient vector of the
 
ambient incident wave, <math>{\bf B}_l</math> is the diffraction transfer
 
matrix of <math>\Delta_l</math> and <math>{\bf T}_{jl}</math> is the coordinate transformation
 
matrix analogous to  (T_elem_deep).
 
 
The calculation of the diffraction transfer matrices is
 
also similar to the infinite depth case. The finite depth
 
Green's function
 
<center><math> (green_d)
 
&G(r,\theta,z;s,\varphi,c)\\ &= \frac{\i}{2} \,
 
\frac{\alpha^2-k^2}{d(\alpha^2-k^2)-\alpha}\, \cosh k(z+d) \cosh k(c+d)
 
\sum_{\nu=-\infty}^{\infty} H_\nu^{(1)}(k r) J_\nu(k s) \mathrm{e}^{\mathrm{i}\nu
 
(\theta - \varphi)}\\ 
 
& \quad + \frac{1}{\pi} \sum_{m=1}^{\infty}
 
\frac{k_m^2+\alpha^2}{d(k_m^2+\alpha^2)-\alpha}\, \cos k_m(z+d) \cos
 
k_m(c+d) \sum_{\nu=-\infty}^{\infty} K_\nu(k_m r) I_\nu(k_m s) \mathrm{e}^{\mathrm{i}\nu
 
(\theta - \varphi)},
 
 
given by [[black75]] and [[fenton78]], needs to be used instead
 
of the infinite depth Green's function  (green_inf).
 
The elements of <math>{\bf B}_j</math> are therefore given by
 
\begin{subequations} (B_elem_d)
 
<center><math>\begin{matrix}
 
({\bf B}_j)_{pq} &= \frac{\i}{2} \,
 
\frac{(\alpha^2-k^2)\cosh^2 kd}{d(\alpha^2-k^2)-\alpha} \int\limits_{\Gamma_j}
 
\cosh k(c+d) J_p(\alpha s) \mathrm{e}^{-\mathrm{i}p \varphi} \varsigma_q^j(\mathbf{\zeta})
 
\mathrm{d}\sigma_\mathbf{\zeta}\\
 
=and= 
 
({\bf B}_j)_{pq} &= \frac{1}{\pi}
 
\frac{(k_m^2+\alpha^2)\cos^2 k_md}{d(k_m^2+\alpha^2)-\alpha}
 
\int\limits_{\Gamma_j} \cos k_m(c+d) I_p(\eta s) \mathrm{e}^{-\mathrm{i}p
 
\varphi} \varsigma_q^j(\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta}
 
\end{matrix}</math></center>
 
\end{subequations}
 
for the propagating and the decaying modes respectively, where
 
<math>\varsigma_q^j(\mathbf{\zeta})</math> is the source strength distribution
 
due to an incident potential of mode <math>q</math> of the form
 
\begin{subequations} (test_modes_d)
 
<center><math>\begin{matrix}
 
\phi_q^{\mathrm{I}}(s,\varphi,c) &=  \frac{\cosh k_m(c+d)}{\cosh kd}
 
H_q^{(1)} (k s) \mathrm{e}^{\mathrm{i}q \varphi}\\
 
=for the propagating modes, and= 
 
\phi_q^{\mathrm{I}}(s,\varphi,c) &= \frac{\cos k_m(c+d)}{\cos kd} K_q
 
(k_m s) \mathrm{e}^{\mathrm{i}q \varphi}
 
\end{matrix}</math></center>
 
\end{subequations}
 
for the decaying modes.
 
 
 
 
 
==Wave forcing of an ice floe of arbitrary geometry==
 
 
The interaction theory which has been developed so far has been
 
for arbitrary bodies. No assumption has been made about the body
 
geometry or its equations of motion. However, we will now use this
 
interaction theory to make calculations for the specific case of ice
 
floes. Ice floes form in vast fields consisting of hundreds if not
 
thousands of individual floes and furthermore  most ice floe fields
 
occur in the deep ocean. For this reason they are ideally suited to
 
the application of the scattering theory we have just developed.
 
Furthermore, the presence of the ice lengthens the wavelength
 
making it more difficult to determine how deep the water must
 
be to be approximately infinite.
 
 
===Mathematical model for ice floes===
 
We will briefly describe the mathematical model which is used to
 
describe ice floes. A more detailed account can be found in
 
[[Squire_review]]. We assume that the ice floe is sufficiently thin
 
that we may apply the shallow draft approximation, which essentially
 
applies the boundary conditions underneath the floe at the water
 
surface. The ice floe is modelled as a thin plate rather than a rigid
 
body since the floe flexure is significant owing to the ice floe
 
geometry.  This model has been applied to a single ice floe by
 
[[JGR02]]. Assuming the ice floe is in contact with the water
 
surface at all times, its displacement
 
<math>W</math> is that of the water surface and <math>W</math> is required to satisfy the linear
 
plate equation in the area occupied by the ice floe <math>\Delta</math>. In analogy to
 
(time), <math>w</math> denotes the time-independent surface displacement
 
(with the same radian frequency as the water velocity potential due to
 
linearity) and the plate equation becomes
 
<center><math> (plate_non)
 
D \, \nabla^4 w - \omega^2 \, \rho_\Delta \, h \, w = \mathrm{i}\, \omega \, \rho
 
\, \phi - \rho \, g \, w, \quad {\bf{x}} \in \Delta,
 
 
with the density of the water <math>\rho</math>, the modulus of rigidity of the
 
ice floe <math>D</math>, its density <math>\rho_\Delta</math> and its
 
thickness <math>h</math>. The right-hand-side of  (plate_non) arises from the
 
linearised Bernoulli equation. It needs to be recalled that
 
<math>\mathbf{x}</math> always denotes a point of the undisturbed water surface.
 
Free edge boundary conditions apply, namely
 
<center><math>
 
\frac{\partial^2 w}{\partial n^2} + \nu \frac{\partial^2 w}{\partial
 
s^2} = 0 \quad =and=  \quad \frac{\partial^3 w}{\partial n^3} + (2 - \nu)
 
\frac{\partial^3 w}{\partial n \partial s^2} = 0, \quad
 
\mathbf{x} \in \partial \Delta,
 
 
where <math>n</math> and <math>s</math> denote the normal and tangential directions on
 
<math>\partial \Delta</math> (where they exist) respectively and <math>\nu</math> is
 
Poisson's ratio.
 
 
Non-dimensional variables (denoted with an overbar) are introduced,
 
<center><math>
 
(\bar{x},\bar{y},\bar{z}) = \frac{1}{a} (x,y,z), \quad \bar{w} =
 
\frac{w}{a}, \quad \bar{\alpha} = a\, \alpha, \quad \bar{\omega} = \omega
 
\sqrt{\frac{a}{g}} \quad =and=  \quad \bar{\phi} = \frac{\phi}{a
 
\sqrt{a g}},
 
 
where <math>a</math> is a length parameter associated with the floe.
 
In non-dimensional variables, the equation for the ice floe
 
(plate_non) reduces to
 
<center><math> (plate_final)
 
\beta \nabla^4 \bar{w} - \bar{\alpha} \gamma \bar{w} = \i
 
\sqrt{\bar{\alpha}}  \bar{\phi} - \bar{w}, \quad
 
\bar{\mathbf{x}} \in \bar{\Delta}
 
 
with
 
<center><math>
 
\beta = \frac{D}{g \rho a^4} \quad =and=  \quad \gamma =
 
\frac{\rho_\Delta h}{ \rho a}.
 
 
The constants <math>\beta</math> and <math>\gamma</math> represent the stiffness and the
 
mass of the plate respectively. For convenience, the overbars will be
 
dropped and non-dimensional variables will be assumed in the sequel.
 
 
The standard boundary-value problem applies to the water.
 
The water velocity potential must satisfy the boundary value problem
 
\begin{subequations} (water)
 
<center><math>\begin{matrix}
 
\nabla^2 \phi &= 0, \; & & \mathbf{y} \in D,\\ 
 
(water_freesurf)
 
\frac{\partial \phi}{\partial z} &= \alpha \phi, \; & &
 
{\bf{x}} \not\in \Delta,\\
 
(water_depth)
 
\sup_{\mathbf{y} \in D} \abs{\phi} &< \infty.
 
\intertext{The linearised kinematic boundary condition is applied under
 
the ice floe,}
 
(water_body)
 
\frac{\partial \phi}{\partial z} &= - \mathrm{i}\sqrt{\alpha} w, \; && {\bf{x}}
 
\in \Delta,
 
\end{matrix}</math></center>
 
and the Sommerfeld radiation condition
 
<center><math>
 
\lim_{\tilde{r} \rightarrow \infty} \sqrt{\tilde{r}} \, \Big(
 
\frac{\partial}{\partial \tilde{r}} - \mathrm{i}k
 
\Big) (\phi - \phi^{\mathrm{In}}) = 0,
 
 
\end{subequations}
 
where <math>\tilde{r}^2=x^2+y^2</math> and <math>k</math> is the wavenumber is imposed.
 
 
Since the numerical convergence will be compared to the finite depth
 
theory later, a formulation for the finite depth problem will be
 
required. However, the differences to the infinite depth
 
formulation are few. For water of constant finite depth <math>d</math>, the volume
 
occupied by the water changes, the vertical dimension being reduced to
 
<math>(-d,0)</math>, (still denoted by <math>D</math>),
 
and the depth condition  (water_depth) is replaced by the bed
 
condition,
 
<center><math>
 
\frac{\partial \phi}{\partial z} = 0, \quad \mathbf{y} \in D,\: z=-d.
 
 
In water of finite depth, the positive real wavenumber <math>k</math> is related
 
to the radian frequency by the dispersion relations  (eq_k).
 
 
 
===The wavelength under the ice floe=== (sec:kappa)
 
For the case of a floating thin plate of shallow draft, which we have
 
used here to model ice floes, waves can propagate under the plate.
 
These
 
waves can be understood by considering an infinite sheet of ice
 
and they satisfy a complex dispersion relation given by
 
[[FoxandSquire]]. In non-dimensional form it states
 
<center><math>
 
\kappa^* \tan \kappa^* d = - \frac{\alpha}{\beta \kappa^{*4} - \gamma
 
\alpha +1},
 
 
where <math>\kappa^*</math> is the wavenumber under the plate. The purely imaginary
 
roots of this dispersion relation correspond to the propagating modes
 
and their absolute value is given as the positive root of
 
<center><math>
 
\kappa \tanh \kappa d = \frac{\alpha}{\beta \kappa^4 - \gamma
 
\alpha + 1}.
 
 
For realistic values of the parameters, the effect of the plate is to
 
make <math>\kappa</math> smaller than <math>k</math> (the open water wavenumber), which
 
increases the wavelength. The effect of the increased wavelength is to
 
increase the depth at which the water may be approximated as
 
infinite.
 
 
===Transformation into an integral equation===
 
The problem for the water is converted to an integral equation in the
 
following way. Let <math>G</math> be the three-dimensional free surface
 
Green's function for water of infinite depth.
 
The Green's function allows the representation of the scattered water
 
velocity potential in the standard way,
 
<center><math> (int_eq)
 
\phi^\mathrm{S}(\mathbf{y}) = \int\limits_{\Gamma}
 
\left( \phi^\mathrm{S} (\mathbf{\zeta}) \, \frac{\partial G}{\partial
 
n_\mathbf{\zeta}} (\mathbf{y};\mathbf{\zeta}) - G
 
(\mathbf{y};\mathbf{\zeta}) \, \frac{\partial
 
\phi^\mathrm{S}}{\partial n_\mathbf{\zeta}} (\mathbf{\zeta}) \right)
 
\mathrm{d}\sigma_\mathbf{\zeta}, \quad \mathbf{y} \in D.
 
 
In the case of a shallow draft, the fact that the Green's function is
 
symmetric and therefore satisfies the free surface boundary condition
 
with respect to the second variable as well can be used to
 
drastically simplify  (int_eq). Due to the linearity of the problem
 
the ambient incident potential can just be added to the equation to obtain the
 
total water velocity potential,
 
<math>\phi=\phi^{\mathrm{I}}+\phi^{\mathrm{S}}</math>. Limiting the result to
 
the water surface leaves the integral equation for the water velocity
 
potential under the ice floe,
 
<center><math> (int_eq_hs)
 
\phi(\mathbf{x}) = \phi^{\mathrm{I}}(\mathbf{x}) +
 
\int\limits_{\Delta} G (\mathbf{x};\mathbf{\xi}) \big( \alpha
 
\phi(\mathbf{\xi}) + \mathrm{i}\sqrt{\alpha} w(\mathbf{\xi}) \big)
 
\mathrm{d}\sigma_\mathbf{\xi}, \quad \mathbf{x} \in \Delta.
 
 
Since the surface displacement of the ice floe appears in this
 
integral equation, it is coupled with the plate equation  (plate_final).
 
A method of solution is discussed in detail by [[JGR02]] but a short
 
outline will be given. The surface displacement of the ice floe is
 
expanded into its modes of vibration by calculating the eigenfunctions
 
and eigenvalues of the biharmonic operator. The integral equation for
 
the potential is then solved for every eigenfunction which gives a
 
corresponding potential to each eigenfunction. The expansion in the
 
eigenfunctions simplifies the biharmonic equation and, by using the
 
orthogonality of the eigenfunctions, a system of equations for the
 
unknown coefficients of the eigenfunction expansion is obtained.
 
 
===The coupled ice floe - water equations===
 
 
Since the operator <math>\nabla^4</math>, subject to the free edge boundary
 
conditions, is self-adjoint a thin plate must possess a set of modes <math>w^k</math>
 
which satisfy the free boundary conditions and the eigenvalue
 
equation
 
<center><math>
 
\nabla^4 w^k = \lambda_k w^k.
 
 
The modes which correspond to different eigenvalues <math>\lambda_k</math> are
 
orthogonal and the eigenvalues are positive and real. While the plate will
 
always have repeated eigenvalues, orthogonal modes can still be found and
 
the modes can be normalised. We therefore assume that the modes are
 
orthonormal, i.e.
 
<center><math>
 
\int\limits_\Delta w^j (\mathbf{\xi}) w^k (\mathbf{\xi})
 
\mathrm{d}\sigma_{\mathbf{\xi}} = \delta _{jk},
 
 
where <math>\delta _{jk}</math> is the Kronecker delta. The eigenvalues <math>\lambda_k</math>
 
have the property that <math>\lambda_k \rightarrow \infty</math> as <math>k \rightarrow
 
\infty</math> and we order the modes by increasing eigenvalue. These modes can be
 
used to expand any function over the wetted surface of the ice floe <math>\Delta</math>.
 
 
We expand the displacement of the floe in a finite number of modes <math>M</math>, i.e.
 
<center><math> (expansion)
 
w(\mathbf{x}) =\sum_{k=1}^{M} c_k w^k (\mathbf{x}).
 
 
>From the linearity of  (int_eq_hs) the potential can be
 
written in the form
 
<center><math> (expansionphi)
 
\phi(\mathbf{x}) =\phi^0(\mathbf{x}) + \sum_{k=1}^{M} c_k \phi^k (\mathbf{x}),
 
 
where <math>\phi^0</math> and <math>\phi^k</math> respectively satisfy the integral equations
 
\begin{subequations} (phi)
 
<center><math> (phi0)
 
\phi^0(\mathbf{x}) = \phi^{\mathrm{I}} (\mathbf{x}) +
 
\int\limits_\Delta \alpha G (\mathbf{x};\mathbf{\xi}) \phi^0
 
(\mathbf{\xi}) d\sigma_\mathbf{\xi}
 
 
and
 
<center><math>  (phii)
 
\phi^k (\mathbf{x}) = \int\limits_{\Delta} G (\mathbf{x};\mathbf{\xi})
 
\left( \alpha \phi^k (\mathbf{\xi}) + \mathrm{i}\sqrt{\alpha} w^k
 
(\mathbf{\xi})\right) \mathrm{d}\sigma_{\mathbf{\xi}}. 
 
 
\end{subequations}
 
The potential <math>\phi^0</math> represents the potential due to the incoming wave
 
assuming that the displacement of the ice floe is zero. The potential
 
<math>\phi^k</math> represents the potential which is generated by the plate
 
vibrating with the <math>k</math>th mode in the absence of any input wave forcing.
 
 
We substitute equations  (expansion) and  (expansionphi) into
 
equation  (plate_final) to obtain
 
<center><math> (expanded)
 
\beta \sum_{k=1}^{M} \lambda_k c_k w^k -\alpha \gamma
 
\sum_{k=1}^{M} c_k w^k = \mathrm{i}\sqrt{\alpha} \big( \phi^0 +
 
\sum_{k=1}^{M} c_k \phi^k \big) - \sum_{k=1}^{M} c_k w^k.
 
 
To solve equation  (expanded) we multiply by <math>w^j</math> and integrate over
 
the plate (i.e. we take the inner product with respect to <math>w^j</math>) taking
 
into account the orthogonality of the modes <math>w^j</math> and obtain
 
<center><math> (final)
 
\beta \lambda_k c_k + \left( 1-\alpha \gamma \right) c_k =
 
\int\limits_{\Delta} \mathrm{i}\sqrt{\alpha} \big( \phi^0 (\mathbf{\xi})
 
+ \sum_{j=1}^{N} c_j \phi^j (\mathbf{\xi}) \big) w^k (\mathbf{\xi})
 
\mathrm{d}\sigma_{\mathbf{\xi}},
 
 
which is a matrix equation in <math>c_k</math>.
 
 
Equation  (final) cannot be solved without determining the modes of
 
vibration of the thin plate <math>w^k</math> (along with the associated
 
eigenvalues <math>\lambda_k</math>) and solving the integral equations
 
(phi). We use the finite element method to
 
determine the modes of vibration \cite[]{Zienkiewicz} and the integral
 
equations  (phi) are solved by a constant panel
 
method \cite[]{Sarp_Isa}. The same set of nodes is used for the finite
 
element method and to define the panels for the integral equation.
 
 
 
===Full diffraction calculation for multiple ice floes===
 
 
The interaction theory is a method to allow more rapid solutions to
 
problems involving multiple bodies. The principle advantage is that the
 
potential is represented in the cylindrical eigenfunctions and
 
therefore fewer unknowns are required. However, every problem which
 
can be solved by the interaction theory can also be solved by applying
 
the full diffraction theory and solving an integral equation over the
 
wetted surface of all the bodies. In this section we will briefly show
 
how this extension can be performed for the ice floe situation. The
 
full diffraction calculation will be used to check the performance and
 
convergence of our interaction theory. Also, because the interaction
 
theory is only valid when the escribed cylinder for each ice floe does
 
not contain any other floe, the full diffraction calculation is
 
required for a very dense arrangement of ice floes.
 
 
We can solve the full diffraction problem for multiple ice floes by the
 
following extension. The displacement of the <math>j</math>th floe is expanded in a
 
finite number of modes <math>M_j</math> (since the number of modes may not
 
necessarily be the same), i.e. 
 
<center><math> (expansion_f)
 
w_j \left( \mathbf{x}\right) =\sum_{k=1}^{M_j} c_{jk} w_j^k (\mathbf{x}). 
 
 
>From the linearity of  (int_eq_hs) the potential can be
 
written in the form
 
<center><math> (expansionphi_f)
 
\phi(\mathbf{x}) =\phi_0(\mathbf{x}) + \sum_{n=1}^{N} \sum_{k=1}^{M_n}
 
c_{nk} \phi_n^k(\mathbf{x}),
 
 
where <math>\phi _{0}</math> and <math>\phi_j^k</math> respectively satisfy the integral equations
 
\begin{subequations} (phi_f)
 
<center><math>  (phi0_f)
 
\phi_j^0 (\mathbf{x}) = \phi^{\mathrm{I}} (\mathbf{x}) + \sum_{n=1}^{N}
 
\int\limits_{\Delta_n} \alpha G (\mathbf{x};\mathbf{\xi})
 
\phi_j^0(\mathbf{\xi}) \mathrm{d}\sigma_{\mathbf{\xi}}
 
 
and
 
<center><math> (phii_f)
 
\phi_j^k(\mathbf{x}) = \sum_{n=1}^{N} \int\limits_{\Delta_n} G
 
(\mathbf{x};\mathbf{\xi}) \left( \alpha \phi_j^k (\mathbf{\xi}) +
 
\i\sqrt{\alpha} w_j^k (\mathbf{\xi})\right) \mathrm{d}\sigma_{\mathbf{\xi}}.
 
 
\end{subequations}
 
The potential <math>\phi_j^{0}</math> represents the potential due the incoming wave
 
assuming that the displacement of the ice floe is zero,
 
<math>\phi_j^k</math> represents the potential which is generated by the <math>j</math>th plate
 
vibrating with the <math>k</math>th mode in the absence of any input wave forcing. It
 
should be noted that <math>\phi_j^k(\mathbf{x})</math> is, in general, non-zero
 
for <math>\mathbf{x}\in \Delta_{n}</math> (since the vibration of the <math>j</math>th
 
plate will result in potential under the <math>n</math>th plate).
 
 
We substitute equations  (expansion_f) and  (expansionphi_f) into
 
equation  (plate_final) to obtain
 
<center><math> (expanded_f)
 
\beta_j \sum_{k=1}^{M_j} \lambda_{jk} c_{jk} w_j^k - \alpha \gamma_j
 
\sum_{k=1}^{M_j} c_{jk} w_j^k = \mathrm{i}\sqrt{\alpha} \,
 
\big( \phi_j^0 + \sum_{n=1}^{N} \sum_{k=1}^{M_n} c_{nk} \phi_n^k \big)
 
- \sum_{k=1}^{M_n} c_{jk} w_j^k. 
 
 
To solve equation  (expanded_f) we multiply by <math>w_j^{l}</math> and integrate
 
over the plate (as before)
 
taking into account the orthogonality of the modes <math>w_j^l</math> and obtain
 
<center><math> (final_f)
 
\beta_j \lambda_{jk} c_{jk} + \big(1- \alpha \gamma_j \big)
 
c_{jk} = \int\limits_{\Delta_j} \mathrm{i}\sqrt{\alpha} \, \big( \phi_0
 
(\mathbf{\xi}) + \sum_{n=1}^{N} \sum_{l=1}^{M_n} c_{nl} \phi_n^l
 
(\mathbf{\xi}) \big) \, w_j^l (\mathbf{\xi}) \mathrm{d}\sigma_{\mathbf{\xi}},
 
 
which holds for all <math>j= 1, \ldots, N</math> and therefore gives a matrix
 
equation for <math>c_{jk}</math>.
 
 
==Numerical Results==
 
 
In this section we will present some calculations using the interaction
 
theory in finite and infinite depth and the full
 
diffraction method in finite and infinite depth.
 
These will be based on calculations for ice floes. We begin with some
 
convergence tests which aim to compare the various methods. It needs
 
to be noted that this comparison is only of numerical nature since the
 
interactions methods as well as the full diffraction calculations
 
are exact in an analytical sense. However, numerical calculations
 
require truncations which affect the different methods in different
 
ways. Especially the dependence on these truncations will be investigated.
 
 
===Convergence Test===
 
We will present some convergence tests that aim to compare the
 
performance of the interaction theory with the full diffraction
 
calculations and to compare the
 
performance of the finite and infinite depth interaction methods in deep water.
 
The comparisons will be conducted for the case of two square ice floes
 
in three different arrangements.
 
In the full diffraction calculation the ice floes
 
are discretised in <math>24 \times 24 = 576</math> elements. For the full diffraction
 
calculation the resulting linear system of equations to be solved is
 
therefore 1152. As will be seen, once the diffraction
 
transfer matrix has been calculated (and saved), the dimension of the
 
linear system of equations to be solved in the interaction method is
 
considerably smaller. It is given by twice the dimension of the
 
diffraction transfer matrix. The most challenging situation for the
 
interaction theory is when the bodies are close together. For this
 
reason we choose the distance such that the escribed circles
 
of the two ice floes just overlap. It must be recalled that the
 
interaction theory is valid as long as the escribed cylinder of a body
 
does not intersect with any other body.
 
 
Both ice floes have non-dimensionalised
 
stiffness <math>\beta = 0.02</math>, mass <math>\gamma = 0.02</math> and Poisson's ratio
 
is chosen as <math>\nu=0.3333</math>. The wavelength of
 
the ambient incident wave is <math>\lambda = 2</math>. Each ice floe has
 
side length 2. The ambient
 
wavefield is of unit amplitude and propagates in the <math>x</math>-direction.
 
Three different arrangements are chosen to compare the results of the
 
finite depth interaction method in deep water and the infinite depth
 
interaction method with the corresponding full diffraction
 
calculations. In the first arrangement the second ice floe is located
 
behind the first, in the second arrangement it is located
 
beside, and the third arrangement it is both
 
beside and behind. The exact positions of the ice floes
 
are given in table (tab:pos).
 
 
\begin{table}
 
\begin{center}
 
\begin{tabular}{@{}ccc@{}}
 
arrangement & <math>O_1</math> & <math>O_2</math>\<center><math>3pt]
 
1 & <math>(-1.4,0)</math> & <math>(1.4,0)</math>\\
 
2 & <math>(0,-1.4)</math> & <math>(0,1.4)</math>\\
 
3 & <math>(-1.4,-0.6)</math> & <math>(1.4,0.6)</math>
 
\end{tabular}
 
\caption{Positions of the ice floes in the different arrangements.} (tab:pos)
 
\end{center}
 
\end{table}
 
 
Figure (fig:tsf) shows the
 
solutions corresponding to the three arrangements in the case of water
 
of infinite depth. To illustrate the effect on the water in the
 
vicinity of the ice floes, the water displacement is also shown.
 
It is interesting
 
to note that the ice floe in front is barely influenced by the
 
floe behind while the motion of the floe behind is quite
 
different from its motion in the absence of the floe in front.
 
 
\begin{figure}
 
\begin{center}
 
 
\includegraphics[width=8.2cm]{int_n11_x0_y28_v5_b02_chi2}\<center><math>0.4cm]
 
\includegraphics[width=8.2cm]{int_n11_x0_y28_v5_b02_chi0}\<center><math>0.4cm]
 
\includegraphics[width=8.2cm]{int_n11_x12_y28_v5_b02_chi2}
 
 
\end{center}
 
\caption{Surface displacement of the ice floes
 
and the water in their vicinity,\newline arrangements 1, 2 and 3.} (fig:tsf)
 
\end{figure}
 
 
To compare the results, a measure of
 
the error from the full diffraction calculation is used. We calculate
 
the full diffraction solution with a sufficient number of points
 
so that we may use it to approximate the exact solution.
 
<center><math>
 
E_2 = \left( \, \int\limits_{\Delta}
 
\big| w_{i}(\mathbf{x})-w_{f}(\mathbf{x}) \big|^2 \mathrm{d}\mathbf{x} \,
 
\right)^{1/2},
 
 
where <math>w_{i}</math> and <math>w_{f}</math> are the solutions of the
 
interaction method and the corresponding full diffraction calculation
 
respectively. It would also be possible to compare other errors, the
 
maximum difference of the solutions for example, but the results are
 
very similar.
 
 
It is worth noting that the finite depth interaction method
 
only converges up to a certain depth if used with the
 
eigenfunction expansion of the finite depth Green's function  (green_d).
 
This is because of the factor
 
<math>\alpha^2-k^2</math> in the term of propagating modes of the Green's
 
function. The Green's function can
 
be rewritten by making use of the dispersion relation  (eq_k)
 
\cite[as suggested by][p. 26, for example]{linton01}
 
and the depth restriction of the finite depth interaction method for
 
bodies of arbitrary geometry can be circumvented.
 
 
The truncation parameters for the interaction methods will
 
now be considered for both finite and infinite depth.
 
The number of propagating modes and angular decaying
 
components are free parameters in both methods. In
 
finite depth, the number of decaying roots of the dispersion relation
 
needs to be chosen while in infinite depth the discretisation of
 
a continuous variable must be selected.
 
In the infinite depth case we are free to choose the number of
 
points as well as the points themselves. In water of finite depth, the depth
 
can also be considered a free parameter as long as it is chosen large
 
enough to account for deep water.
 
 
Truncating the infinite sums in the eigenfunction expansion of the
 
outgoing water velocity potential for infinite depth with
 
truncation parameters <math>T_H</math> and <math>T_K</math> and discretising the integration
 
by defining a set of nodes, <math>0\leq\eta_1 < \ldots < \eta_m < \ldots <
 
\eta_{_{T_R}}<math>, with weights </math>h_m</math>, the potential for infinite depth
 
can be approximated by
 
<center><math>
 
\phi (r,\theta,z) &=  \mathrm{e}^{\alpha z} \sum_{\nu = -
 
T_H}^{T_H} A_{0\nu} H_\nu^{(1)} (\alpha r) \mathrm{e}^{\mathrm{i}\nu
 
\theta}\\
 
&\quad + \sum_{m=1}^{T_R} h_m \, \psi (z,\eta_m) \sum_{\nu = -
 
T_K}^{T_K} A_{\nu} (\eta_m) K_\nu (\eta_m r) \mathrm{e}^{\mathrm{i}\nu \theta}.
 
 
In the following, the integration weights are chosen to be </math>h_m =
 
1/2\,(\eta_{m+1}-\eta_{m-1})<math>, </math>m=2, \ldots, T_R-1<math> and </math>h_1 =
 
\eta_2-\eta_1<math> as well as </math>h_{_{T_R}} =
 
\eta_{_{T_R}}-\eta_{_{T_R-1}}<math>, which corresponds to the mid-point
 
quadrature rule.
 
Different quadrature rules such as Gaussian quadrature
 
could be considered. Although in general this would lead to better
 
results, the mid-point rule allows a clever
 
choice of the discretisation points so that the convergence with
 
Gaussian quadrature is no better.
 
In finite depth, the analogous truncation leads to
 
<center><math>
 
\phi (r,\theta,z) &= \frac{\cosh k (z+d)}{\cosh kd} \sum_{\nu = -
 
T_H}^{T_H} A_{0\nu} H_\nu^{(1)} (k r) \mathrm{e}^{\mathrm{i}\nu \theta}\\
 
& \quad + \sum_{m = 1}^{T_R} \frac{\cos k_m(z+d)}{\cos k_m d}
 
\sum_{\nu = - T_K}^{T_K} A_{m\nu} K_\nu (k_m r) \mathrm{e}^{\mathrm{i}\nu \theta}.
 
 
In both cases, the dimension of the diffraction transfer matrix,
 
<math>\mathbf{B}</math>, is given by <math>2 \, T_H+1+T_R \, (2 \, T_K+1)</math>.
 
 
Since the choice of the number of propagating
 
modes and angular decaying components affects the finite and
 
infinite depth methods in similar ways, the dependence on these
 
parameters will not be further presented. Thorough convergence tests
 
have shown that in the settings investigated here, it is sufficient to
 
choose <math>T_H</math> to be 11 and <math>T_K</math> to be 5. Further increasing these
 
parameter values does not result in smaller errors (as compared
 
to the full diffraction calculation with 576 elements per floe).
 
We will now compare the convergence of the infinite depth and
 
the finite depth methods if <math>T_H</math> and <math>T_K</math> are
 
fixed (with the previously mentioned values) and <math>T_R</math> is varied. To be able to
 
compare the results, the discretisation of the continuous variable
 
will always be the same for fixed <math>T_R</math> and these are
 
shown in table (tab:discr).
 
It should be noted that if only one node is used the integration
 
weight is chosen to be 1.
 
 
\begin{table}
 
\begin{center}
 
\begin{tabular}{@{}cl@{}}
 
<math>T_R</math> & discretisation of <math>\eta</math>\<center><math>3pt]
 
1 & \{ 2.1 \}\\
 
2 & \{ 1.2, 2.7 \}\\
 
3 & \{ 0.8, 1.8, 3.0 \}\\
 
4 & \{ 0.4, 1.4, 2.2, 3.2 \}\\
 
5 & \{ 0.2, 1.0, 1.8, 3.0, 4.6 \}
 
\end{tabular}
 
\caption{The different discretisations used in the convergence tests.} (tab:discr)
 
\end{center}
 
\end{table}
 
 
Figures (fig:behind), (fig:beside) and (fig:shifted) show the convergence for arrangement 1, arrangement 2 and arrangement 3, respectively, for
 
the infinite depth method and the finite depth method with depth 2
 
(plot (a)) and depth 4 (plot (b)).
 
Since the ice floes are located beside each other
 
in arrangement 2 the average errors are the same for both floes.
 
As can be seen from figures (fig:behind), (fig:beside) and
 
(fig:shifted) the convergence of the infinite depth method
 
is similar to that of the finite depth method. Used with depth 2, the
 
convergence of the finite depth method is generally better than that
 
of the infinite depth method while used with depth 4, the infinite depth
 
method achieves the better results. Tests with other depths show that
 
the performance of the finite depth method decreases with increasing
 
water depth as expected. In general, since the wavelength is 2, a depth
 
of <math>d=2</math> should approximate infinite depth and hence there is no
 
advantage to using the infinite depth theory. However, as mentioned
 
previously, for certain situations such as ice floes it is not necessarily
 
true that <math>d=2</math> will approximate infinite depth. 
 
 
\begin{figure}
 
\begin{center}
 
\begin{tabular}{p{0.46\columnwidth}p{0.02\columnwidth}p{0.46\columnwidth}} 
 
\includegraphics[height=0.38\columnwidth]{behind_d2}&&
 
\includegraphics[height=0.38\columnwidth]{behind_d4}
 
\end{tabular}
 
\caption{Development of the errors as <math>T_R</math> is increased in
 
arrangement 1.} (fig:behind)
 
\end{center}
 
\end{figure}
 
 
\begin{figure}
 
\begin{center}
 
\begin{tabular}{p{0.46\columnwidth}p{0.02\columnwidth}p{0.46\columnwidth}} 
 
\includegraphics[height=0.38\columnwidth]{beside_d2}&&
 
\includegraphics[height=0.38\columnwidth]{beside_d4}
 
\end{tabular}
 
\caption{Development of the errors as <math>T_R</math> is increased in
 
arrangement 2.} (fig:beside)
 
\end{center}
 
\end{figure}
 
 
\begin{figure}
 
\begin{center}
 
\begin{tabular}{p{0.46\columnwidth}p{0.02\columnwidth}p{0.46\columnwidth}} 
 
\includegraphics[height=0.38\columnwidth]{shifted_d2}&&
 
\includegraphics[height=0.38\columnwidth]{shifted_d4}
 
\end{tabular}
 
\caption{Development of the errors as <math>T_R</math> is increased in
 
arrangement 3.} (fig:shifted)
 
\end{center}
 
\end{figure}
 
 
===Multiple ice floe results===
 
We will now present results for multiple ice floes of different
 
geometries and in different arrangements on water of infinite depth.
 
We choose the floe arrangements arbitrarily, since there are
 
no known special ice floe arrangements, such as those that give
 
rise to resonances in the infinite limit.
 
In all plots, the wavelength <math>\lambda</math> has been chosen to
 
be <math>2</math>, the stiffness <math>\beta</math> and the mass <math>\gamma</math> of the ice
 
floes to be 0.02 and Poisson's ratio <math>\nu</math> is <math>0.3333</math>. The ambient
 
wavefield of amplitude 1 propagates in
 
the positive direction of the <math>x</math>-axis, thus it travels from left to
 
right in the plots. 
 
 
Figure (fig:int_arb) shows the
 
displacements of multiple interacting ice floes of different shapes and
 
in different arrangements. Since square elements have been used to
 
represent the floes, non-rectangular geometries are approximated.
 
All ice floes have an area of 4 and the escribing circles do not
 
intersect with any of the other ice floes.
 
The plots show the displacement of the ice floes at time <math>t=0</math>.
 
 
\begin{figure}
 
\begin{center}
 
\begin{tabular}{p{0.46\columnwidth}p{0.03\columnwidth}p{0.46\columnwidth}}
 
\includegraphics[width=0.45\columnwidth]{mult_n15_tr_tr_tr_tr_in} &&
 
\includegraphics[width=0.45\columnwidth]{mult_n15_tr_tr_tr_tr_out}\<center><math>0.2cm]
 
\includegraphics[width=0.45\columnwidth]{mult_n15_rh_rh_rh} &&
 
\includegraphics[width=0.45\columnwidth]{mult_n15_sq_sq_sq_sq_sq}\\
 
\end{tabular}
 
\end{center}
 
\caption{Surface displacement of interacting ice floes of different geometries.} (fig:int_arb)
 
\end{figure}
 
 
 
 
 
==Summary==
 
The finite depth interaction theory developed by
 
[[kagemoto86]] has been extended to water of infinite
 
depth. Furthermore, using the eigenfunction
 
expansion of the infinite depth free surface Green's function we have
 
been able to calculate the diffraction transfer matrices for bodies of
 
arbitrary geometry. We also showed how the diffraction transfer
 
matrices can be calculated efficiently for different orientations of
 
the body.
 
 
The convergence of the infinite depth interaction method is similar to
 
that of the finite depth method. Generally, it can be said that the
 
greater the water depth in the finite depth method the poorer its
 
performance. Since bodies in the water can change the water depth
 
which is required to allow the water to be approximated as infinitely
 
deep (ice floes for example) it is recommendable to use the infinite
 
depth method if the water depth may be considered
 
infinite. Furthermore, the infinite depth method requires the infinite
 
depth single diffraction solutions which are easier to
 
compute than the finite depth solutions.
 
It is also possible that the
 
convergence of the infinite depth method may be further improved
 
by a novel to optimisation of the discretisation of the continuous variable.
 
 
 
 
 
 
 
 
We extend the finite depth interaction theory of [[kagemoto86]] to
 
water of infinite depth and bodies of arbitrary geometry. The sum
 
over the discrete roots of the dispersion equation in the finite depth
 
theory becomes
 
an integral in the infinite depth theory. This means that the infinite
 
dimensional diffraction
 
transfer matrix
 
in the finite depth theory must be replaced by an integral
 
operator. In the numerical solution of the equations, this
 
integral operator is approximated by a sum and a linear system
 
of equations is obtained. We also show how the calculations
 
of the diffraction transfer matrix for bodies of arbitrary
 
geometry developed by [[goo90]] can be extended to
 
infinite depth, and how the diffraction transfer matrix for rotated bodies can
 
be easily calculated. This interaction theory is applied to the wave forcing
 
of multiple ice floes and a method to solve
 
the full diffraction problem in this case is presented. Convergence
 
studies comparing the interaction method with the full diffraction
 
calculations and the finite and infinite depth interaction methods are
 
carried out. 
 
 
 
 
==Introduction==
 
The scattering of water waves by floating or submerged
 
bodies is of wide practical importance.
 
Although the problem is non-linear, if the
 
wave amplitude is sufficiently small, the
 
problem can be linearised. The linear problem is still the basis of
 
the engineering design of most off-shore structures and is the standard
 
model of geophysical phenomena such as the wave forcing of ice floes.
 
While analytic solutions have been found for simplified problems
 
(especially for simple geometries or in two dimensions) the full
 
three-dimensional linear diffraction problem can only be solved by
 
numerical methods involving the discretisation of the body's surface.
 
The resulting linear system of equations has a dimension equal to the
 
number of unknowns used in the discretisation of the body.
 
 
If more than one body is present, all bodies scatter the
 
incoming waves. Therefore, the scattered wave from one body is
 
incident upon all the others and, given that they are not too far apart,
 
this notably changes the total incident wave upon them.
 
Therefore, the diffraction calculation must be conducted
 
for all bodies simultaneously. Since each body must be discretised this
 
can lead to a very large number of unknowns. However, the scattered
 
wavefield can be represented in an eigenfunction basis with a
 
comparatively small number of unknowns. If we can express the problem in this
 
basis, using what is known as an {\em Interaction Theory},
 
the number of unknowns can be much reduced, especially if there is a large
 
number of bodies.
 
 
The first interaction theory that was not based on an approximation was
 
the interaction theory
 
developed by [[kagemoto86]]. Kagemoto and Yue found an exact
 
algebraic method to solve the linear wave problem for vertically
 
non-overlapping bodies in water of finite depth.  The only restriction
 
of their theory was that the smallest escribed circle for each body must not
 
overlap any other body. The interaction of the bodies was
 
accounted for by taking the scattered wave of each body to be the
 
incident wave upon all other bodies (in addition to the ambient
 
incident wave). Furthermore, since the cylindrical eigenfunction expansions
 
are local, these were mapped from one body to another using
 
Graf's addition theorem for Bessel functions.
 
Doing this for all bodies,
 
\citeauthor{kagemoto86} were able to solve for the
 
coefficients of the scattered wavefields of all bodies simultaneously.
 
The only difficulty with this method was that the
 
solutions of the single diffraction problems had to be available in
 
the cylindrical eigenfunction expansion of an outgoing
 
wave. \citeauthor{kagemoto86} therefore only solved for axisymmetric
 
bodies because the single diffraction solution for axisymmetric
 
bodies was
 
available in the required representation. 
 
 
The extension of the Kagemoto and Yue scattering theory to bodies of
 
arbitrary geometry was performed by [[goo90]] who found a general
 
way to solve the single diffraction problem in the required
 
cylindrical eigenfunction representation. They used a
 
representation of the finite depth free surface Green's function in
 
the eigenfunction expansion of cylindrical outgoing waves
 
centred at an arbitrary point of the water surface (above the
 
body's mean centre position in this case). This Green's function was
 
presented by
 
[[black75]] and further investigated by [[fenton78]]
 
who corrected some statements about the Green's function which Black had
 
made. This Green's function is
 
based on the cylindrical
 
eigenfunction expansion of the finite depth free surface Green's
 
function given by [[john2]]. The results of
 
\citeauthor{goo90} were recently used by [[chakrabarti00]] to
 
solve for arrays of cylinders which can be divided into modules.
 
 
The development of the Kagemoto and Yue interaction theory was
 
motivated by problems
 
in off-shore engineering. However, the theory can also be applied
 
to the geophysical problem of wave scattering by ice floes. At the
 
interface of the open and frozen ocean an interfacial region known
 
as the Marginal Ice Zone (MIZ) forms. The MIZ largely controls
 
the interaction of the open and frozen ocean, especially the interaction
 
through wave processors. This is because
 
the MIZ consists of vast fields
 
of ice floes whose size is comparable to the dominant wavelength, which
 
means that it strongly scatters incoming waves. A method of
 
solving for the wave response of a single ice floe of arbitrary
 
geometry in water of infinite depth was presented by
 
[[JGR02]]. The ice floe was modelled as a floating, flexible
 
thin plate and its motion was expanded in the free plate modes of
 
vibration. Converting the problem for the water into an integral
 
equation and substituting the free modes, a system of equations for
 
the coefficients in the modal expansion was obtained. However,
 
to understand wave propagation and scattering in the MIZ we need to
 
understand the way in which large numbers of interacting
 
ice floes scatter waves. For this reason, we require an interaction
 
theory. While the Kagemoto and Yue interaction theory could be used,
 
their theory requires that the water depth is finite.
 
While the water depth in the Marginal Ice Zone varies,
 
it is generally located far from shore above the deep ocean. This means
 
that the finite depth must be chosen large in order to be able to
 
apply their theory.  Furthermore, when ocean waves propagate beneath
 
an ice floe the wavelength is increased so that it becomes more
 
difficult to make the
 
water depth sufficiently deep that it may be approximated as infinite. For this
 
reason, in this paper we develop  the equivalent interaction theory to
 
Kagemoto and Yue's in infinite depth. Also, because of
 
the complicated geometry of an ice floe, this interaction theory is
 
for bodies of arbitrary geometry. 
 
 
In the first part of this paper Kagemoto and Yue's interaction theory
 
is extended to water of infinite depth. We represent the incident and
 
scattered potentials in the cylindrical eigenfunction expansions
 
and we use an analogous infinite depth Green's function
 
to the one used by \citeauthor{goo90} \cite[given by][]{malte03}. We
 
show how the infinite
 
depth diffraction transfer matrices can be obtained with the use of
 
this Green's function and we illustrate how the rotation of a body
 
about its mean centre position in the plane can be accounted for without
 
recalculating the diffraction transfer matrix.
 
 
In the second part of the paper, using \citeauthor{JGR02}'s single
 
floe result,
 
the full diffraction calculation for the motion and scattering from many
 
interacting ice floes is calculated and presented. For two square
 
interacting ice floes the convergence of the method obtained from the
 
developed interaction theory is compared to the result of the full
 
diffraction calculation. The solutions of more than two interacting
 
ice floes and of other shapes in different arrangements are presented as well.
 
We also compare the convergence of the finite depth and infinite
 
depth methods in
 
deep water.
 
  
 
\section{The extension of Kagemoto and Yue's interaction
 
\section{The extension of Kagemoto and Yue's interaction
Line 2,207: Line 587:
 
easily transformed into a classical <math>(N \, n)</math>-dimensional linear system of
 
easily transformed into a classical <math>(N \, n)</math>-dimensional linear system of
 
equations.
 
equations.
 
 
==Finite Depth Interaction Theory==
 
 
We will compare the performance of the infinite depth interaction theory
 
with the equivalent theory for finite
 
depth. As we have stated previously, the finite depth theory was
 
developed by [[kagemoto86]] and extended to bodies of arbitrary
 
geometry by [[goo90]]. We will briefly present this theory in
 
our notation and the comparisons will be made in a later section.
 
 
In water of constant finite depth <math>d</math>, the scattered potential of a body
 
<math>\Delta_j</math> can be expanded in cylindrical eigenfunctions,
 
<center><math> (basisrep_out_d)
 
\phi_j^\mathrm{S} (r_j,\theta_j,z) &= \frac{\cosh k(z+d)}{\cosh kd}
 
\sum_{\nu = - \infty}^{\infty} A_{0 \nu}^j H_\nu^{(1)} (k r_j) \mathrm{e}^{\mathrm{i}\nu
 
\theta_j}\\
 
&\quad + \sum_{m=1}^{\infty} \frac{\cos k_m (z+d)}{\cos kd} \sum_{\nu = -
 
\infty}^{\infty} A_{m \nu}^j K_\nu (k_m r_j) \mathrm{e}^{\mathrm{i}\nu
 
\theta_j},
 
</math></center>
 
with discrete coefficients <math>A_{m \nu}^j</math>. The positive wavenumber <math>k</math>
 
is related to <math>\alpha</math> by the dispersion relation
 
<center><math> (eq_k)
 
\alpha = k \tanh k d,
 
</math></center>
 
and the values of <math>k_m</math>, <math>m>0</math>, are given as positive real roots of
 
the dispersion relation
 
<center><math> (eq_k_m)
 
\alpha + k_m \tan k_m d = 0.
 
</math></center>
 
The incident potential upon body <math>\Delta_j</math> can be also be expanded in
 
cylindrical eigenfunctions,
 
<center><math> (basisrep_in_d)
 
\phi_j^\mathrm{I} (r_j,\theta_j,z) &= \frac{\cosh k(z+d)}{\cosh kd}
 
\sum_{\mu = - \infty}^{\infty} D_{0 \mu}^j J_\mu (k r_j) \mathrm{e}^{\mathrm{i}\mu
 
\theta_j}\\
 
& \quad + \sum_{m=1}^{\infty} \frac{\cos k_m (z+d)}{\cos kd} \sum_{\mu = -
 
\infty}^{\infty} D_{m\mu}^j I_\mu (k_m r_j) \mathrm{e}^{\mathrm{i}\mu
 
\theta_j},
 
</math></center>
 
with discrete coefficients <math>D_{m\mu}^j</math>. A system of equations for the
 
coefficients of the scattered wavefields for the bodies are derived
 
in an analogous way to the infinite depth case. The derivation is
 
simpler because all the coefficients are discrete and the
 
diffraction transfer operator can be represented by an
 
infinite dimensional matrix.
 
Truncating the infinite dimensional matrix as well as the
 
coefficient vectors appropriately, the resulting system of
 
equations is given by 
 
<center><math>
 
{\bf a}_l = {\bf B}_l \Big( {\bf d}_l^\mathrm{In} +
 
\sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \trans {\bf T}_{jl} \,
 
{\bf a}_j \Big), \quad  l=1, \ldots, N,
 
</math></center>
 
where <math>{\bf a}_l</math> is the coefficient vector of the scattered
 
wave, <math>{\bf d}_l^\mathrm{In}</math> is the coefficient vector of the
 
ambient incident wave, <math>{\bf B}_l</math> is the diffraction transfer
 
matrix of <math>\Delta_l</math> and <math>{\bf T}_{jl}</math> is the coordinate transformation
 
matrix analogous to  (T_elem_deep).
 
 
The calculation of the diffraction transfer matrices is
 
also similar to the infinite depth case. The finite depth
 
Green's function
 
<center><math> (green_d)
 
&G(r,\theta,z;s,\varphi,c)\\ &= \frac{\i}{2} \,
 
\frac{\alpha^2-k^2}{d(\alpha^2-k^2)-\alpha}\, \cosh k(z+d) \cosh k(c+d)
 
\sum_{\nu=-\infty}^{\infty} H_\nu^{(1)}(k r) J_\nu(k s) \mathrm{e}^{\mathrm{i}\nu
 
(\theta - \varphi)}\\ 
 
& \quad + \frac{1}{\pi} \sum_{m=1}^{\infty}
 
\frac{k_m^2+\alpha^2}{d(k_m^2+\alpha^2)-\alpha}\, \cos k_m(z+d) \cos
 
k_m(c+d) \sum_{\nu=-\infty}^{\infty} K_\nu(k_m r) I_\nu(k_m s) \mathrm{e}^{\mathrm{i}\nu
 
(\theta - \varphi)},
 
</math></center>
 
given by [[black75]] and [[fenton78]], needs to be used instead
 
of the infinite depth Green's function  (green_inf).
 
The elements of <math>{\bf B}_j</math> are therefore given by
 
\begin{subequations} (B_elem_d)
 
<center><math>\begin{matrix}
 
({\bf B}_j)_{pq} &= \frac{\i}{2} \,
 
\frac{(\alpha^2-k^2)\cosh^2 kd}{d(\alpha^2-k^2)-\alpha} \int\limits_{\Gamma_j}
 
\cosh k(c+d) J_p(\alpha s) \mathrm{e}^{-\mathrm{i}p \varphi} \varsigma_q^j(\mathbf{\zeta})
 
\mathrm{d}\sigma_\mathbf{\zeta}\\
 
=and= 
 
({\bf B}_j)_{pq} &= \frac{1}{\pi}
 
\frac{(k_m^2+\alpha^2)\cos^2 k_md}{d(k_m^2+\alpha^2)-\alpha}
 
\int\limits_{\Gamma_j} \cos k_m(c+d) I_p(\eta s) \mathrm{e}^{-\mathrm{i}p
 
\varphi} \varsigma_q^j(\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta}
 
\end{matrix}</math></center>
 
\end{subequations}
 
for the propagating and the decaying modes respectively, where
 
<math>\varsigma_q^j(\mathbf{\zeta})</math> is the source strength distribution
 
due to an incident potential of mode <math>q</math> of the form
 
\begin{subequations} (test_modes_d)
 
<center><math>\begin{matrix}
 
\phi_q^{\mathrm{I}}(s,\varphi,c) &=  \frac{\cosh k_m(c+d)}{\cosh kd}
 
H_q^{(1)} (k s) \mathrm{e}^{\mathrm{i}q \varphi}\\
 
=for the propagating modes, and= 
 
\phi_q^{\mathrm{I}}(s,\varphi,c) &= \frac{\cos k_m(c+d)}{\cos kd} K_q
 
(k_m s) \mathrm{e}^{\mathrm{i}q \varphi}
 
\end{matrix}</math></center>
 
\end{subequations}
 
for the decaying modes.
 
 
 
 
 
==Wave forcing of an ice floe of arbitrary geometry==
 
 
The interaction theory which has been developed so far has been
 
for arbitrary bodies. No assumption has been made about the body
 
geometry or its equations of motion. However, we will now use this
 
interaction theory to make calculations for the specific case of ice
 
floes. Ice floes form in vast fields consisting of hundreds if not
 
thousands of individual floes and furthermore  most ice floe fields
 
occur in the deep ocean. For this reason they are ideally suited to
 
the application of the scattering theory we have just developed.
 
Furthermore, the presence of the ice lengthens the wavelength
 
making it more difficult to determine how deep the water must
 
be to be approximately infinite.
 
 
===Mathematical model for ice floes===
 
We will briefly describe the mathematical model which is used to
 
describe ice floes. A more detailed account can be found in
 
[[Squire_review]]. We assume that the ice floe is sufficiently thin
 
that we may apply the shallow draft approximation, which essentially
 
applies the boundary conditions underneath the floe at the water
 
surface. The ice floe is modelled as a thin plate rather than a rigid
 
body since the floe flexure is significant owing to the ice floe
 
geometry.  This model has been applied to a single ice floe by
 
[[JGR02]]. Assuming the ice floe is in contact with the water
 
surface at all times, its displacement
 
<math>W</math> is that of the water surface and <math>W</math> is required to satisfy the linear
 
plate equation in the area occupied by the ice floe <math>\Delta</math>. In analogy to
 
(time), <math>w</math> denotes the time-independent surface displacement
 
(with the same radian frequency as the water velocity potential due to
 
linearity) and the plate equation becomes
 
<center><math> (plate_non)
 
D \, \nabla^4 w - \omega^2 \, \rho_\Delta \, h \, w = \mathrm{i}\, \omega \, \rho
 
\, \phi - \rho \, g \, w, \quad {\bf{x}} \in \Delta,
 
</math></center>
 
with the density of the water <math>\rho</math>, the modulus of rigidity of the
 
ice floe <math>D</math>, its density <math>\rho_\Delta</math> and its
 
thickness <math>h</math>. The right-hand-side of  (plate_non) arises from the
 
linearised Bernoulli equation. It needs to be recalled that
 
<math>\mathbf{x}</math> always denotes a point of the undisturbed water surface.
 
Free edge boundary conditions apply, namely
 
<center><math>
 
\frac{\partial^2 w}{\partial n^2} + \nu \frac{\partial^2 w}{\partial
 
s^2} = 0 \quad =and=  \quad \frac{\partial^3 w}{\partial n^3} + (2 - \nu)
 
\frac{\partial^3 w}{\partial n \partial s^2} = 0, \quad
 
\mathbf{x} \in \partial \Delta,
 
</math></center>
 
where <math>n</math> and <math>s</math> denote the normal and tangential directions on
 
<math>\partial \Delta</math> (where they exist) respectively and <math>\nu</math> is
 
Poisson's ratio.
 
 
Non-dimensional variables (denoted with an overbar) are introduced,
 
<center><math>
 
(\bar{x},\bar{y},\bar{z}) = \frac{1}{a} (x,y,z), \quad \bar{w} =
 
\frac{w}{a}, \quad \bar{\alpha} = a\, \alpha, \quad \bar{\omega} = \omega
 
\sqrt{\frac{a}{g}} \quad =and=  \quad \bar{\phi} = \frac{\phi}{a
 
\sqrt{a g}},
 
</math></center>
 
where <math>a</math> is a length parameter associated with the floe.
 
In non-dimensional variables, the equation for the ice floe
 
(plate_non) reduces to
 
<center><math> (plate_final)
 
\beta \nabla^4 \bar{w} - \bar{\alpha} \gamma \bar{w} = \i
 
\sqrt{\bar{\alpha}}  \bar{\phi} - \bar{w}, \quad
 
\bar{\mathbf{x}} \in \bar{\Delta}
 
</math></center>
 
with
 
<center><math>
 
\beta = \frac{D}{g \rho a^4} \quad =and=  \quad \gamma =
 
\frac{\rho_\Delta h}{ \rho a}.
 
</math></center>
 
The constants <math>\beta</math> and <math>\gamma</math> represent the stiffness and the
 
mass of the plate respectively. For convenience, the overbars will be
 
dropped and non-dimensional variables will be assumed in the sequel.
 
 
The standard boundary-value problem applies to the water.
 
The water velocity potential must satisfy the boundary value problem
 
\begin{subequations} (water)
 
<center><math>\begin{matrix}
 
\nabla^2 \phi &= 0, \; & & \mathbf{y} \in D,\\ 
 
(water_freesurf)
 
\frac{\partial \phi}{\partial z} &= \alpha \phi, \; & &
 
{\bf{x}} \not\in \Delta,\\
 
(water_depth)
 
\sup_{\mathbf{y} \in D} \abs{\phi} &< \infty.
 
\intertext{The linearised kinematic boundary condition is applied under
 
the ice floe,}
 
(water_body)
 
\frac{\partial \phi}{\partial z} &= - \mathrm{i}\sqrt{\alpha} w, \; && {\bf{x}}
 
\in \Delta,
 
\end{matrix}</math></center>
 
and the Sommerfeld radiation condition
 
<center><math>
 
\lim_{\tilde{r} \rightarrow \infty} \sqrt{\tilde{r}} \, \Big(
 
\frac{\partial}{\partial \tilde{r}} - \mathrm{i}k
 
\Big) (\phi - \phi^{\mathrm{In}}) = 0,
 
</math></center>
 
\end{subequations}
 
where <math>\tilde{r}^2=x^2+y^2</math> and <math>k</math> is the wavenumber is imposed.
 
 
Since the numerical convergence will be compared to the finite depth
 
theory later, a formulation for the finite depth problem will be
 
required. However, the differences to the infinite depth
 
formulation are few. For water of constant finite depth <math>d</math>, the volume
 
occupied by the water changes, the vertical dimension being reduced to
 
<math>(-d,0)</math>, (still denoted by <math>D</math>),
 
and the depth condition  (water_depth) is replaced by the bed
 
condition,
 
<center><math>
 
\frac{\partial \phi}{\partial z} = 0, \quad \mathbf{y} \in D,\: z=-d.
 
</math></center>
 
In water of finite depth, the positive real wavenumber <math>k</math> is related
 
to the radian frequency by the dispersion relations  (eq_k).
 
 
 
===The wavelength under the ice floe=== (sec:kappa)
 
For the case of a floating thin plate of shallow draft, which we have
 
used here to model ice floes, waves can propagate under the plate.
 
These
 
waves can be understood by considering an infinite sheet of ice
 
and they satisfy a complex dispersion relation given by
 
[[FoxandSquire]]. In non-dimensional form it states
 
<center><math>
 
\kappa^* \tan \kappa^* d = - \frac{\alpha}{\beta \kappa^{*4} - \gamma
 
\alpha +1},
 
</math></center>
 
where <math>\kappa^*</math> is the wavenumber under the plate. The purely imaginary
 
roots of this dispersion relation correspond to the propagating modes
 
and their absolute value is given as the positive root of
 
<center><math>
 
\kappa \tanh \kappa d = \frac{\alpha}{\beta \kappa^4 - \gamma
 
\alpha + 1}.
 
</math></center>
 
For realistic values of the parameters, the effect of the plate is to
 
make <math>\kappa</math> smaller than <math>k</math> (the open water wavenumber), which
 
increases the wavelength. The effect of the increased wavelength is to
 
increase the depth at which the water may be approximated as
 
infinite.
 
 
===Transformation into an integral equation===
 
The problem for the water is converted to an integral equation in the
 
following way. Let <math>G</math> be the three-dimensional free surface
 
Green's function for water of infinite depth.
 
The Green's function allows the representation of the scattered water
 
velocity potential in the standard way,
 
<center><math> (int_eq)
 
\phi^\mathrm{S}(\mathbf{y}) = \int\limits_{\Gamma}
 
\left( \phi^\mathrm{S} (\mathbf{\zeta}) \, \frac{\partial G}{\partial
 
n_\mathbf{\zeta}} (\mathbf{y};\mathbf{\zeta}) - G
 
(\mathbf{y};\mathbf{\zeta}) \, \frac{\partial
 
\phi^\mathrm{S}}{\partial n_\mathbf{\zeta}} (\mathbf{\zeta}) \right)
 
\mathrm{d}\sigma_\mathbf{\zeta}, \quad \mathbf{y} \in D.
 
</math></center>
 
In the case of a shallow draft, the fact that the Green's function is
 
symmetric and therefore satisfies the free surface boundary condition
 
with respect to the second variable as well can be used to
 
drastically simplify  (int_eq). Due to the linearity of the problem
 
the ambient incident potential can just be added to the equation to obtain the
 
total water velocity potential,
 
<math>\phi=\phi^{\mathrm{I}}+\phi^{\mathrm{S}}</math>. Limiting the result to
 
the water surface leaves the integral equation for the water velocity
 
potential under the ice floe,
 
<center><math> (int_eq_hs)
 
\phi(\mathbf{x}) = \phi^{\mathrm{I}}(\mathbf{x}) +
 
\int\limits_{\Delta} G (\mathbf{x};\mathbf{\xi}) \big( \alpha
 
\phi(\mathbf{\xi}) + \mathrm{i}\sqrt{\alpha} w(\mathbf{\xi}) \big)
 
\mathrm{d}\sigma_\mathbf{\xi}, \quad \mathbf{x} \in \Delta.
 
</math></center>
 
Since the surface displacement of the ice floe appears in this
 
integral equation, it is coupled with the plate equation  (plate_final).
 
A method of solution is discussed in detail by [[JGR02]] but a short
 
outline will be given. The surface displacement of the ice floe is
 
expanded into its modes of vibration by calculating the eigenfunctions
 
and eigenvalues of the biharmonic operator. The integral equation for
 
the potential is then solved for every eigenfunction which gives a
 
corresponding potential to each eigenfunction. The expansion in the
 
eigenfunctions simplifies the biharmonic equation and, by using the
 
orthogonality of the eigenfunctions, a system of equations for the
 
unknown coefficients of the eigenfunction expansion is obtained.
 
 
===The coupled ice floe - water equations===
 
 
Since the operator <math>\nabla^4</math>, subject to the free edge boundary
 
conditions, is self-adjoint a thin plate must possess a set of modes <math>w^k</math>
 
which satisfy the free boundary conditions and the eigenvalue
 
equation
 
<center><math>
 
\nabla^4 w^k = \lambda_k w^k.
 
</math></center>
 
The modes which correspond to different eigenvalues <math>\lambda_k</math> are
 
orthogonal and the eigenvalues are positive and real. While the plate will
 
always have repeated eigenvalues, orthogonal modes can still be found and
 
the modes can be normalised. We therefore assume that the modes are
 
orthonormal, i.e.
 
<center><math>
 
\int\limits_\Delta w^j (\mathbf{\xi}) w^k (\mathbf{\xi})
 
\mathrm{d}\sigma_{\mathbf{\xi}} = \delta _{jk},
 
</math></center>
 
where <math>\delta _{jk}</math> is the Kronecker delta. The eigenvalues <math>\lambda_k</math>
 
have the property that <math>\lambda_k \rightarrow \infty</math> as </math>k \rightarrow
 
\infty<math> and we order the modes by increasing eigenvalue. These modes can be
 
used to expand any function over the wetted surface of the ice floe <math>\Delta</math>.
 
 
We expand the displacement of the floe in a finite number of modes <math>M</math>, i.e.
 
<center><math> (expansion)
 
w(\mathbf{x}) =\sum_{k=1}^{M} c_k w^k (\mathbf{x}).
 
</math></center>
 
>From the linearity of  (int_eq_hs) the potential can be
 
written in the form
 
<center><math> (expansionphi)
 
\phi(\mathbf{x}) =\phi^0(\mathbf{x}) + \sum_{k=1}^{M} c_k \phi^k (\mathbf{x}),
 
</math></center>
 
where <math>\phi^0</math> and <math>\phi^k</math> respectively satisfy the integral equations
 
\begin{subequations} (phi)
 
<center><math> (phi0)
 
\phi^0(\mathbf{x}) = \phi^{\mathrm{I}} (\mathbf{x}) +
 
\int\limits_\Delta \alpha G (\mathbf{x};\mathbf{\xi}) \phi^0
 
(\mathbf{\xi}) d\sigma_\mathbf{\xi}
 
</math></center>
 
and
 
<center><math>  (phii)
 
\phi^k (\mathbf{x}) = \int\limits_{\Delta} G (\mathbf{x};\mathbf{\xi})
 
\left( \alpha \phi^k (\mathbf{\xi}) + \mathrm{i}\sqrt{\alpha} w^k
 
(\mathbf{\xi})\right) \mathrm{d}\sigma_{\mathbf{\xi}}. 
 
</math></center>
 
\end{subequations}
 
The potential <math>\phi^0</math> represents the potential due to the incoming wave
 
assuming that the displacement of the ice floe is zero. The potential
 
<math>\phi^k</math> represents the potential which is generated by the plate
 
vibrating with the <math>k</math>th mode in the absence of any input wave forcing.
 
 
We substitute equations  (expansion) and  (expansionphi) into
 
equation  (plate_final) to obtain
 
<center><math> (expanded)
 
\beta \sum_{k=1}^{M} \lambda_k c_k w^k -\alpha \gamma
 
\sum_{k=1}^{M} c_k w^k = \mathrm{i}\sqrt{\alpha} \big( \phi^0 +
 
\sum_{k=1}^{M} c_k \phi^k \big) - \sum_{k=1}^{M} c_k w^k.
 
</math></center>
 
To solve equation  (expanded) we multiply by <math>w^j</math> and integrate over
 
the plate (i.e. we take the inner product with respect to <math>w^j</math>) taking
 
into account the orthogonality of the modes <math>w^j</math> and obtain
 
<center><math> (final)
 
\beta \lambda_k c_k + \left( 1-\alpha \gamma \right) c_k =
 
\int\limits_{\Delta} \mathrm{i}\sqrt{\alpha} \big( \phi^0 (\mathbf{\xi})
 
+ \sum_{j=1}^{N} c_j \phi^j (\mathbf{\xi}) \big) w^k (\mathbf{\xi})
 
\mathrm{d}\sigma_{\mathbf{\xi}},
 
</math></center>
 
which is a matrix equation in <math>c_k</math>.
 
 
Equation  (final) cannot be solved without determining the modes of
 
vibration of the thin plate <math>w^k</math> (along with the associated
 
eigenvalues <math>\lambda_k</math>) and solving the integral equations
 
(phi). We use the finite element method to
 
determine the modes of vibration \cite[]{Zienkiewicz} and the integral
 
equations  (phi) are solved by a constant panel
 
method \cite[]{Sarp_Isa}. The same set of nodes is used for the finite
 
element method and to define the panels for the integral equation.
 
 
 
===Full diffraction calculation for multiple ice floes===
 
 
The interaction theory is a method to allow more rapid solutions to
 
problems involving multiple bodies. The principle advantage is that the
 
potential is represented in the cylindrical eigenfunctions and
 
therefore fewer unknowns are required. However, every problem which
 
can be solved by the interaction theory can also be solved by applying
 
the full diffraction theory and solving an integral equation over the
 
wetted surface of all the bodies. In this section we will briefly show
 
how this extension can be performed for the ice floe situation. The
 
full diffraction calculation will be used to check the performance and
 
convergence of our interaction theory. Also, because the interaction
 
theory is only valid when the escribed cylinder for each ice floe does
 
not contain any other floe, the full diffraction calculation is
 
required for a very dense arrangement of ice floes.
 
 
We can solve the full diffraction problem for multiple ice floes by the
 
following extension. The displacement of the <math>j</math>th floe is expanded in a
 
finite number of modes <math>M_j</math> (since the number of modes may not
 
necessarily be the same), i.e. 
 
<center><math> (expansion_f)
 
w_j \left( \mathbf{x}\right) =\sum_{k=1}^{M_j} c_{jk} w_j^k (\mathbf{x}). 
 
</math></center>
 
>From the linearity of  (int_eq_hs) the potential can be
 
written in the form
 
<center><math> (expansionphi_f)
 
\phi(\mathbf{x}) =\phi_0(\mathbf{x}) + \sum_{n=1}^{N} \sum_{k=1}^{M_n}
 
c_{nk} \phi_n^k(\mathbf{x}),
 
</math></center>
 
where <math>\phi _{0}</math> and <math>\phi_j^k</math> respectively satisfy the integral equations
 
\begin{subequations} (phi_f)
 
<center><math>  (phi0_f)
 
\phi_j^0 (\mathbf{x}) = \phi^{\mathrm{I}} (\mathbf{x}) + \sum_{n=1}^{N}
 
\int\limits_{\Delta_n} \alpha G (\mathbf{x};\mathbf{\xi})
 
\phi_j^0(\mathbf{\xi}) \mathrm{d}\sigma_{\mathbf{\xi}}
 
</math></center>
 
and
 
<center><math> (phii_f)
 
\phi_j^k(\mathbf{x}) = \sum_{n=1}^{N} \int\limits_{\Delta_n} G
 
(\mathbf{x};\mathbf{\xi}) \left( \alpha \phi_j^k (\mathbf{\xi}) +
 
\i\sqrt{\alpha} w_j^k (\mathbf{\xi})\right) \mathrm{d}\sigma_{\mathbf{\xi}}.
 
</math></center>
 
\end{subequations}
 
The potential <math>\phi_j^{0}</math> represents the potential due the incoming wave
 
assuming that the displacement of the ice floe is zero,
 
<math>\phi_j^k</math> represents the potential which is generated by the <math>j</math>th plate
 
vibrating with the <math>k</math>th mode in the absence of any input wave forcing. It
 
should be noted that <math>\phi_j^k(\mathbf{x})</math> is, in general, non-zero
 
for <math>\mathbf{x}\in \Delta_{n}</math> (since the vibration of the <math>j</math>th
 
plate will result in potential under the <math>n</math>th plate).
 
 
We substitute equations  (expansion_f) and  (expansionphi_f) into
 
equation  (plate_final) to obtain
 
<center><math> (expanded_f)
 
\beta_j \sum_{k=1}^{M_j} \lambda_{jk} c_{jk} w_j^k - \alpha \gamma_j
 
\sum_{k=1}^{M_j} c_{jk} w_j^k = \mathrm{i}\sqrt{\alpha} \,
 
\big( \phi_j^0 + \sum_{n=1}^{N} \sum_{k=1}^{M_n} c_{nk} \phi_n^k \big)
 
- \sum_{k=1}^{M_n} c_{jk} w_j^k. 
 
</math></center>
 
To solve equation  (expanded_f) we multiply by <math>w_j^{l}</math> and integrate
 
over the plate (as before)
 
taking into account the orthogonality of the modes <math>w_j^l</math> and obtain
 
<center><math> (final_f)
 
\beta_j \lambda_{jk} c_{jk} + \big(1- \alpha \gamma_j \big)
 
c_{jk} = \int\limits_{\Delta_j} \mathrm{i}\sqrt{\alpha} \, \big( \phi_0
 
(\mathbf{\xi}) + \sum_{n=1}^{N} \sum_{l=1}^{M_n} c_{nl} \phi_n^l
 
(\mathbf{\xi}) \big) \, w_j^l (\mathbf{\xi}) \mathrm{d}\sigma_{\mathbf{\xi}},
 
</math></center>
 
which holds for all <math>j= 1, \ldots, N</math> and therefore gives a matrix
 
equation for <math>c_{jk}</math>.
 
 
==Numerical Results==
 
 
In this section we will present some calculations using the interaction
 
theory in finite and infinite depth and the full
 
diffraction method in finite and infinite depth.
 
These will be based on calculations for ice floes. We begin with some
 
convergence tests which aim to compare the various methods. It needs
 
to be noted that this comparison is only of numerical nature since the
 
interactions methods as well as the full diffraction calculations
 
are exact in an analytical sense. However, numerical calculations
 
require truncations which affect the different methods in different
 
ways. Especially the dependence on these truncations will be investigated.
 
 
===Convergence Test===
 
We will present some convergence tests that aim to compare the
 
performance of the interaction theory with the full diffraction
 
calculations and to compare the
 
performance of the finite and infinite depth interaction methods in deep water.
 
The comparisons will be conducted for the case of two square ice floes
 
in three different arrangements.
 
In the full diffraction calculation the ice floes
 
are discretised in <math>24 \times 24 = 576</math> elements. For the full diffraction
 
calculation the resulting linear system of equations to be solved is
 
therefore 1152. As will be seen, once the diffraction
 
transfer matrix has been calculated (and saved), the dimension of the
 
linear system of equations to be solved in the interaction method is
 
considerably smaller. It is given by twice the dimension of the
 
diffraction transfer matrix. The most challenging situation for the
 
interaction theory is when the bodies are close together. For this
 
reason we choose the distance such that the escribed circles
 
of the two ice floes just overlap. It must be recalled that the
 
interaction theory is valid as long as the escribed cylinder of a body
 
does not intersect with any other body.
 
 
Both ice floes have non-dimensionalised
 
stiffness <math>\beta = 0.02</math>, mass <math>\gamma = 0.02</math> and Poisson's ratio
 
is chosen as <math>\nu=0.3333</math>. The wavelength of
 
the ambient incident wave is <math>\lambda = 2</math>. Each ice floe has
 
side length 2. The ambient
 
wavefield is of unit amplitude and propagates in the <math>x</math>-direction.
 
Three different arrangements are chosen to compare the results of the
 
finite depth interaction method in deep water and the infinite depth
 
interaction method with the corresponding full diffraction
 
calculations. In the first arrangement the second ice floe is located
 
behind the first, in the second arrangement it is located
 
beside, and the third arrangement it is both
 
beside and behind. The exact positions of the ice floes
 
are given in table (tab:pos).
 
 
\begin{table}
 
\begin{center}
 
\begin{tabular}{@{}ccc@{}}
 
arrangement & <math>O_1</math> & <math>O_2</math>\<center><math>3pt]
 
1 & <math>(-1.4,0)</math> & <math>(1.4,0)</math>\\
 
2 & <math>(0,-1.4)</math> & <math>(0,1.4)</math>\\
 
3 & <math>(-1.4,-0.6)</math> & <math>(1.4,0.6)</math>
 
\end{tabular}
 
\caption{Positions of the ice floes in the different arrangements.} (tab:pos)
 
\end{center}
 
\end{table}
 
 
Figure (fig:tsf) shows the
 
solutions corresponding to the three arrangements in the case of water
 
of infinite depth. To illustrate the effect on the water in the
 
vicinity of the ice floes, the water displacement is also shown.
 
It is interesting
 
to note that the ice floe in front is barely influenced by the
 
floe behind while the motion of the floe behind is quite
 
different from its motion in the absence of the floe in front.
 
 
\begin{figure}
 
\begin{center}
 
 
\includegraphics[width=8.2cm]{int_n11_x0_y28_v5_b02_chi2}\<center><math>0.4cm]
 
\includegraphics[width=8.2cm]{int_n11_x0_y28_v5_b02_chi0}\<center><math>0.4cm]
 
\includegraphics[width=8.2cm]{int_n11_x12_y28_v5_b02_chi2}
 
 
\end{center}
 
\caption{Surface displacement of the ice floes
 
and the water in their vicinity,\newline arrangements 1, 2 and 3.} (fig:tsf)
 
\end{figure}
 
 
To compare the results, a measure of
 
the error from the full diffraction calculation is used. We calculate
 
the full diffraction solution with a sufficient number of points
 
so that we may use it to approximate the exact solution.
 
<center><math>
 
E_2 = \left( \, \int\limits_{\Delta}
 
\big| w_{i}(\mathbf{x})-w_{f}(\mathbf{x}) \big|^2 \mathrm{d}\mathbf{x} \,
 
\right)^{1/2},
 
</math></center>
 
where <math>w_{i}</math> and <math>w_{f}</math> are the solutions of the
 
interaction method and the corresponding full diffraction calculation
 
respectively. It would also be possible to compare other errors, the
 
maximum difference of the solutions for example, but the results are
 
very similar.
 
 
It is worth noting that the finite depth interaction method
 
only converges up to a certain depth if used with the
 
eigenfunction expansion of the finite depth Green's function  (green_d).
 
This is because of the factor
 
<math>\alpha^2-k^2</math> in the term of propagating modes of the Green's
 
function. The Green's function can
 
be rewritten by making use of the dispersion relation  (eq_k)
 
\cite[as suggested by][p. 26, for example]{linton01}
 
and the depth restriction of the finite depth interaction method for
 
bodies of arbitrary geometry can be circumvented.
 
 
The truncation parameters for the interaction methods will
 
now be considered for both finite and infinite depth.
 
The number of propagating modes and angular decaying
 
components are free parameters in both methods. In
 
finite depth, the number of decaying roots of the dispersion relation
 
needs to be chosen while in infinite depth the discretisation of
 
a continuous variable must be selected.
 
In the infinite depth case we are free to choose the number of
 
points as well as the points themselves. In water of finite depth, the depth
 
can also be considered a free parameter as long as it is chosen large
 
enough to account for deep water.
 
 
Truncating the infinite sums in the eigenfunction expansion of the
 
outgoing water velocity potential for infinite depth with
 
truncation parameters <math>T_H</math> and <math>T_K</math> and discretising the integration
 
by defining a set of nodes, <math>0\leq\eta_1 < \ldots < \eta_m < \ldots <
 
\eta_{_{T_R}}<math>, with weights </math>h_m</math>, the potential for infinite depth
 
can be approximated by
 
<center><math>
 
\phi (r,\theta,z) &=  \mathrm{e}^{\alpha z} \sum_{\nu = -
 
T_H}^{T_H} A_{0\nu} H_\nu^{(1)} (\alpha r) \mathrm{e}^{\mathrm{i}\nu
 
\theta}\\
 
&\quad + \sum_{m=1}^{T_R} h_m \, \psi (z,\eta_m) \sum_{\nu = -
 
T_K}^{T_K} A_{\nu} (\eta_m) K_\nu (\eta_m r) \mathrm{e}^{\mathrm{i}\nu \theta}.
 
</math></center>
 
In the following, the integration weights are chosen to be <math>h_m =
 
1/2\,(\eta_{m+1}-\eta_{m-1})<math>, </math>m=2, \ldots, T_R-1<math> and </math>h_1 =
 
\eta_2-\eta_1<math> as well as </math>h_{_{T_R}} =
 
\eta_{_{T_R}}-\eta_{_{T_R-1}}</math>, which corresponds to the mid-point
 
quadrature rule.
 
Different quadrature rules such as Gaussian quadrature
 
could be considered. Although in general this would lead to better
 
results, the mid-point rule allows a clever
 
choice of the discretisation points so that the convergence with
 
Gaussian quadrature is no better.
 
In finite depth, the analogous truncation leads to
 
<center><math>
 
\phi (r,\theta,z) &= \frac{\cosh k (z+d)}{\cosh kd} \sum_{\nu = -
 
T_H}^{T_H} A_{0\nu} H_\nu^{(1)} (k r) \mathrm{e}^{\mathrm{i}\nu \theta}\\
 
& \quad + \sum_{m = 1}^{T_R} \frac{\cos k_m(z+d)}{\cos k_m d}
 
\sum_{\nu = - T_K}^{T_K} A_{m\nu} K_\nu (k_m r) \mathrm{e}^{\mathrm{i}\nu \theta}.
 
</math></center>
 
In both cases, the dimension of the diffraction transfer matrix,
 
<math>\mathbf{B}</math>, is given by <math>2 \, T_H+1+T_R \, (2 \, T_K+1)</math>.
 
 
Since the choice of the number of propagating
 
modes and angular decaying components affects the finite and
 
infinite depth methods in similar ways, the dependence on these
 
parameters will not be further presented. Thorough convergence tests
 
have shown that in the settings investigated here, it is sufficient to
 
choose <math>T_H</math> to be 11 and <math>T_K</math> to be 5. Further increasing these
 
parameter values does not result in smaller errors (as compared
 
to the full diffraction calculation with 576 elements per floe).
 
We will now compare the convergence of the infinite depth and
 
the finite depth methods if <math>T_H</math> and <math>T_K</math> are
 
fixed (with the previously mentioned values) and <math>T_R</math> is varied. To be able to
 
compare the results, the discretisation of the continuous variable
 
will always be the same for fixed <math>T_R</math> and these are
 
shown in table (tab:discr).
 
It should be noted that if only one node is used the integration
 
weight is chosen to be 1.
 
 
\begin{table}
 
\begin{center}
 
\begin{tabular}{@{}cl@{}}
 
<math>T_R</math> & discretisation of <math>\eta</math>\<center><math>3pt]
 
1 & \{ 2.1 \}\\
 
2 & \{ 1.2, 2.7 \}\\
 
3 & \{ 0.8, 1.8, 3.0 \}\\
 
4 & \{ 0.4, 1.4, 2.2, 3.2 \}\\
 
5 & \{ 0.2, 1.0, 1.8, 3.0, 4.6 \}
 
\end{tabular}
 
\caption{The different discretisations used in the convergence tests.} (tab:discr)
 
\end{center}
 
\end{table}
 
 
Figures (fig:behind), (fig:beside) and (fig:shifted) show the convergence for arrangement 1, arrangement 2 and arrangement 3, respectively, for
 
the infinite depth method and the finite depth method with depth 2
 
(plot (a)) and depth 4 (plot (b)).
 
Since the ice floes are located beside each other
 
in arrangement 2 the average errors are the same for both floes.
 
As can be seen from figures (fig:behind), (fig:beside) and
 
(fig:shifted) the convergence of the infinite depth method
 
is similar to that of the finite depth method. Used with depth 2, the
 
convergence of the finite depth method is generally better than that
 
of the infinite depth method while used with depth 4, the infinite depth
 
method achieves the better results. Tests with other depths show that
 
the performance of the finite depth method decreases with increasing
 
water depth as expected. In general, since the wavelength is 2, a depth
 
of <math>d=2</math> should approximate infinite depth and hence there is no
 
advantage to using the infinite depth theory. However, as mentioned
 
previously, for certain situations such as ice floes it is not necessarily
 
true that <math>d=2</math> will approximate infinite depth. 
 
 
\begin{figure}
 
\begin{center}
 
\begin{tabular}{p{0.46\columnwidth}p{0.02\columnwidth}p{0.46\columnwidth}} 
 
\includegraphics[height=0.38\columnwidth]{behind_d2}&&
 
\includegraphics[height=0.38\columnwidth]{behind_d4}
 
\end{tabular}
 
\caption{Development of the errors as <math>T_R</math> is increased in
 
arrangement 1.} (fig:behind)
 
\end{center}
 
\end{figure}
 
 
\begin{figure}
 
\begin{center}
 
\begin{tabular}{p{0.46\columnwidth}p{0.02\columnwidth}p{0.46\columnwidth}} 
 
\includegraphics[height=0.38\columnwidth]{beside_d2}&&
 
\includegraphics[height=0.38\columnwidth]{beside_d4}
 
\end{tabular}
 
\caption{Development of the errors as <math>T_R</math> is increased in
 
arrangement 2.} (fig:beside)
 
\end{center}
 
\end{figure}
 
 
\begin{figure}
 
\begin{center}
 
\begin{tabular}{p{0.46\columnwidth}p{0.02\columnwidth}p{0.46\columnwidth}} 
 
\includegraphics[height=0.38\columnwidth]{shifted_d2}&&
 
\includegraphics[height=0.38\columnwidth]{shifted_d4}
 
\end{tabular}
 
\caption{Development of the errors as <math>T_R</math> is increased in
 
arrangement 3.} (fig:shifted)
 
\end{center}
 
\end{figure}
 
 
===Multiple ice floe results===
 
We will now present results for multiple ice floes of different
 
geometries and in different arrangements on water of infinite depth.
 
We choose the floe arrangements arbitrarily, since there are
 
no known special ice floe arrangements, such as those that give
 
rise to resonances in the infinite limit.
 
In all plots, the wavelength <math>\lambda</math> has been chosen to
 
be <math>2</math>, the stiffness <math>\beta</math> and the mass <math>\gamma</math> of the ice
 
floes to be 0.02 and Poisson's ratio <math>\nu</math> is <math>0.3333</math>. The ambient
 
wavefield of amplitude 1 propagates in
 
the positive direction of the <math>x</math>-axis, thus it travels from left to
 
right in the plots. 
 
 
Figure (fig:int_arb) shows the
 
displacements of multiple interacting ice floes of different shapes and
 
in different arrangements. Since square elements have been used to
 
represent the floes, non-rectangular geometries are approximated.
 
All ice floes have an area of 4 and the escribing circles do not
 
intersect with any of the other ice floes.
 
The plots show the displacement of the ice floes at time <math>t=0</math>.
 
 
\begin{figure}
 
\begin{center}
 
\begin{tabular}{p{0.46\columnwidth}p{0.03\columnwidth}p{0.46\columnwidth}}
 
\includegraphics[width=0.45\columnwidth]{mult_n15_tr_tr_tr_tr_in} &&
 
\includegraphics[width=0.45\columnwidth]{mult_n15_tr_tr_tr_tr_out}\<center><math>0.2cm]
 
\includegraphics[width=0.45\columnwidth]{mult_n15_rh_rh_rh} &&
 
\includegraphics[width=0.45\columnwidth]{mult_n15_sq_sq_sq_sq_sq}\\
 
\end{tabular}
 
\end{center}
 
\caption{Surface displacement of interacting ice floes of different geometries.} (fig:int_arb)
 
\end{figure}
 
 
 
 
 
==Summary==
 
The finite depth interaction theory developed by
 
[[kagemoto86]] has been extended to water of infinite
 
depth. Furthermore, using the eigenfunction
 
expansion of the infinite depth free surface Green's function we have
 
been able to calculate the diffraction transfer matrices for bodies of
 
arbitrary geometry. We also showed how the diffraction transfer
 
matrices can be calculated efficiently for different orientations of
 
the body.
 
 
The convergence of the infinite depth interaction method is similar to
 
that of the finite depth method. Generally, it can be said that the
 
greater the water depth in the finite depth method the poorer its
 
performance. Since bodies in the water can change the water depth
 
which is required to allow the water to be approximated as infinitely
 
deep (ice floes for example) it is recommendable to use the infinite
 
depth method if the water depth may be considered
 
infinite. Furthermore, the infinite depth method requires the infinite
 
depth single diffraction solutions which are easier to
 
compute than the finite depth solutions.
 
It is also possible that the
 
convergence of the infinite depth method may be further improved
 
by a novel to optimisation of the discretisation of the continuous variable.</math>
 

Revision as of 02:07, 7 June 2006

This is an interaction theory which provides the exact solution (i.e. it is not based on a Wide Spacing Approximation). The theory uses the Cylindrical Eigenfunction Expansion and Graf's Addition Theorem to represent the potential in local coordinates. The incident and scattered potential of each body are then related by the associated Diffraction Transfer Matrix.

The basic idea is as follows: The scattered potential of each body is represented in the Cylindrical Eigenfunction Expansion associated with the local coordinates centred at the mean centre position of the body. Using Graf's Addition Theorem, the scattered potential of all bodies (given in their local coordinates) can be mapped to an incident potential associated with the coordiates of all other bodies. Doing this, the incident potential of each body (which is given by the ambient incident potential plus the scattered potentials of all other bodies) is given in the Cylindrical Eigenfunction Expansion associated with its local coordinates. Using the Diffraction Transfer Matrix, which relates the incident and scattered potential of each body in isolation, a system of equations for the coefficients of the scattered potentials of all bodies is obtained.

The theory is described in Kagemoto and Yue 1986 and in Peter and Meylan 2004.

\section{The extension of Kagemoto and Yue's interaction theory to bodies of arbitrary shape in water of infinite depth}

kagemoto86 developed an interaction theory for vertically non-overlapping axisymmetric structures in water of finite depth. While their theory was valid for bodies of arbitrary geometry, they did not develop all the necessary details to apply the theory to arbitrary bodies. The only requirements to apply this scattering theory is that the bodies are vertically non-overlapping and that the smallest cylinder which completely contains each body does not intersect with any other body. In this section we will extend their theory to bodies of arbitrary geometry in water of infinite depth. The extension of \citeauthor{kagemoto86}'s finite depth interaction theory to bodies of arbitrary geometry was accomplished by goo90.


The interaction theory begins by representing the scattered potential of each body in the cylindrical eigenfunction expansion. Furthermore, the incoming potential is also represented in the cylindrical eigenfunction expansion. The operator which maps the incoming and outgoing representation is called the diffraction transfer matrix and is different for each body. Since these representations are local to each body, a mapping of the eigenfunction representations between different bodies is required. This operator is called the coordinate transformation matrix.

The cylindrical eigenfunction expansions will be introduced before we derive a system of equations for the coefficients of the scattered wavefields. Analogously to kagemoto86, we represent the scattered wavefield of each body as an incoming wave upon all other bodies. The addition of the ambient incident wave yields the complete incident potential and with the use of diffraction transfer matrices which relate the coefficients of the incident potential to those of the scattered wavefield a system of equations for the unknown coefficients of the scattered wavefields of all bodies is derived.


Eigenfunction expansion of the potential

The equations of motion for the water are derived from the linearised inviscid theory. Under the assumption of irrotational motion the velocity vector field of the water can be written as the gradient field of a scalar velocity potential [math]\displaystyle{ \Phi }[/math]. Assuming that the motion is time-harmonic with the radian frequency [math]\displaystyle{ \omega }[/math] the velocity potential can be expressed as the real part of a complex quantity,

[math]\displaystyle{ (time) \Phi(\mathbf{y},t) = \Re \{\phi (\mathbf{y}) \mathrm{e}^{-\mathrm{i}\omega t} \}. }[/math]

To simplify notation, [math]\displaystyle{ \mathbf{y} = (x,y,z) }[/math] will always denote a point in the water, which is assumed infinitely deep, while [math]\displaystyle{ \mathbf{x} }[/math] will always denote a point of the undisturbed water surface assumed at [math]\displaystyle{ z=0 }[/math].

The problem consists of [math]\displaystyle{ N }[/math] vertically non-overlapping bodies, denoted by [math]\displaystyle{ \Delta_j }[/math], which are sufficiently far apart that there is no intersection of the smallest cylinder which contains each body with any other body. Each body is subject to an incident wavefield which is incoming, responds to this wavefield and produces a scattered wave field which is outgoing. Both the incident and scattered potential corresponding to these wavefields can be represented in the cylindrical eigenfunction expansion valid outside of the escribed cylinder of the body. Let [math]\displaystyle{ (r_j,\theta_j,z) }[/math] be the local cylindrical coordinates of the [math]\displaystyle{ j }[/math]th body, [math]\displaystyle{ \Delta_j }[/math], [math]\displaystyle{ j \in \{1, \ldots , N\} }[/math], and </math>\alpha = \omega^2/g[math]\displaystyle{ where }[/math]g[math]\displaystyle{ is the acceleration due to gravity. Figure (fig:floe_tri) shows these coordinate systems for two bodies. \begin{figure} \begin{center} \includegraphics[height=5.5cm]{floe_tri} \caption{Plan view of the relation between two bodies.} (fig:floe_tri) \end{center} \end{figure} The scattered potential of body \lt math\gt \Delta_j }[/math] can be expanded in cylindrical eigenfunctions,

[math]\displaystyle{ (basisrep_out) \phi_j^\mathrm{S} (r_j,\theta_j,z) &= \mathrm{e}^{\alpha z} \sum_{\nu = - \infty}^{\infty} A_{0 \nu}^j H_\nu^{(1)} (\alpha r_j) \mathrm{e}^{\mathrm{i}\nu \theta_j}\\ &\quad + \int\limits_0^{\infty} \big( \cos \eta z + \frac{\alpha}{\eta} \sin \eta z \big) \sum_{\nu = - \infty}^{\infty} A_{\nu}^j (\eta) K_\nu (\eta r_j) \mathrm{e}^{\mathrm{i}\nu \theta_j} \mathrm{d}\eta, }[/math]

where the coefficients [math]\displaystyle{ A_{0 \nu}^j }[/math] for the propagating modes are discrete and the coefficients [math]\displaystyle{ A_{\nu}^j (\cdot) }[/math] for the decaying modes are functions. [math]\displaystyle{ H_\nu^{(1)} }[/math] and [math]\displaystyle{ K_\nu }[/math] are the Hankel function of the first kind and the modified Bessel function of the second kind respectively, both of order [math]\displaystyle{ \nu }[/math] as defined in abr_ste. The incident potential upon body [math]\displaystyle{ \Delta_j }[/math] can be expanded in cylindrical eigenfunctions,

[math]\displaystyle{ (basisrep_in) \phi_j^\mathrm{I} (r_j,\theta_j,z) &= \mathrm{e}^{\alpha z} \sum_{\mu = - \infty}^{\infty} D_{0 \mu}^j J_\mu (\alpha r_j) \mathrm{e}^{\mathrm{i}\mu \theta_j}\\ & \quad + \int\limits_0^{\infty} \big( \cos \eta z + \frac{\alpha}{\eta} \sin \eta z \big) \sum_{\mu = - \infty}^{\infty} D_{\mu}^j (\eta) I_\mu (\eta r_j) \mathrm{e}^{\mathrm{i}\mu \theta_j} \mathrm{d}\eta, }[/math]

where the coefficients [math]\displaystyle{ D_{0 \mu}^j }[/math] for the propagating modes are discrete and the coefficients [math]\displaystyle{ D_{\mu}^j (\cdot) }[/math] for the decaying modes are functions. [math]\displaystyle{ J_\mu }[/math] and [math]\displaystyle{ I_\mu }[/math] are the Bessel function and the modified Bessel function respectively, both of the first kind and order [math]\displaystyle{ \mu }[/math]. To simplify the notation, from now on [math]\displaystyle{ \psi(z,\eta) }[/math] will denote the vertical eigenfunctions corresponding to the decaying modes,

[math]\displaystyle{ \psi(z,\eta) = \cos \eta z + \alpha / \eta \, \sin \eta z. }[/math]

The interaction in water of infinite depth

Following the ideas of kagemoto86, a system of equations for the unknown coefficients and coefficient functions of the scattered wavefields will be developed. This system of equations is based on transforming the scattered potential of [math]\displaystyle{ \Delta_j }[/math] into an incident potential upon [math]\displaystyle{ \Delta_l }[/math] ([math]\displaystyle{ j \neq l }[/math]). Doing this for all bodies simultaneously, and relating the incident and scattered potential for each body, a system of equations for the unknown coefficients will be developed.

The scattered potential [math]\displaystyle{ \phi_j^{\mathrm{S}} }[/math] of body [math]\displaystyle{ \Delta_j }[/math] needs to be represented in terms of the incident potential [math]\displaystyle{ \phi_l^{\mathrm{I}} }[/math] upon [math]\displaystyle{ \Delta_l }[/math], [math]\displaystyle{ j \neq l }[/math]. From figure (fig:floe_tri) we can see that this can be accomplished by using Graf's addition theorem for Bessel functions given in \citet[eq. 9.1.79]{abr_ste}, \begin{subequations} (transf)

[math]\displaystyle{ \begin{matrix} (transf_h) H_\nu^{(1)}(\alpha r_j) \mathrm{e}^{\mathrm{i}\nu (\theta_j - \vartheta_{jl})} &= \sum_{\mu = - \infty}^{\infty} H^{(1)}_{\nu + \mu} (\alpha R_{jl}) \, J_\mu (\alpha r_l) \mathrm{e}^{\mathrm{i}\mu (\pi - \theta_l + \vartheta_{jl})}, \quad j \neq l,\\ (transf_k) K_\nu(\eta r_j) \mathrm{e}^{\mathrm{i}\nu (\theta_j - \vartheta_{jl})} &= \sum_{\mu = - \infty}^{\infty} K_{\nu + \mu} (\eta R_{jl}) \, I_\mu (\eta r_l) \mathrm{e}^{\mathrm{i}\mu (\pi - \theta_l + \vartheta_{jl})}, \quad j \neq l, \end{matrix} }[/math]

\end{subequations} which is valid provided that [math]\displaystyle{ r_l \lt R_{jl} }[/math]. This limitation only requires that the escribed cylinder of each body [math]\displaystyle{ \Delta_l }[/math] does not enclose any other origin [math]\displaystyle{ O_j }[/math] ([math]\displaystyle{ j \neq l }[/math]). However, the expansion of the scattered and incident potential in cylindrical eigenfunctions is only valid outside the escribed cylinder of each body. Therefore the condition that the escribed cylinder of each body [math]\displaystyle{ \Delta_l }[/math] does not enclose any other origin [math]\displaystyle{ O_j }[/math] ([math]\displaystyle{ j \neq l }[/math]) is superseded by the more rigorous restriction that the escribed cylinder of each body may not contain any other body. Making use of the equations (transf) the scattered potential of [math]\displaystyle{ \Delta_j }[/math] can be expressed in terms of the incident potential upon [math]\displaystyle{ \Delta_l }[/math],

[math]\displaystyle{ \begin{matrix} \phi_j^{\mathrm{S}} (r_l,\theta_l,z) &= \mathrm{e}^{\alpha z} \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j \sum_{\mu = -\infty}^{\infty} H_{\nu - \mu}^{(1)} (\alpha R_{jl}) J_\mu (\alpha r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l} \mathrm{e}^{\mathrm{i}(\nu-\mu) \vartheta_{jl}}\\ & \quad + \int\limits_0^{\infty} \psi (z,\eta) \sum_{\nu = - \infty}^{\infty} A_{\nu}^j (\eta) \sum_{\mu = -\infty}^{\infty} (-1)^\mu K_{\nu-\mu} (\eta R_{jl}) I_\mu (\eta r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l} \mathrm{e}^{\mathrm{i}(\nu-\mu) \vartheta_{jl}} \mathrm{d}\eta\\ &= \mathrm{e}^{\alpha z} \sum_{\mu = -\infty}^{\infty} \Big[ \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j H_{\nu - \mu}^{(1)} (\alpha R_{jl}) \mathrm{e}^{\i (\nu-\mu) \vartheta_{jl}} \Big] J_\mu (\alpha r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l}\\ & + \int\limits_0^{\infty} \psi (z,\eta) \sum_{\mu = -\infty}^{\infty} \Big[ \sum_{\nu = - \infty}^{\infty} A_{\nu}^j (\eta) (-1)^\mu K_{\nu-\mu} (\eta R_{jl}) \mathrm{e}^{\mathrm{i}(\nu-\mu) \vartheta_{jl}} \Big] I_\mu (\eta r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l} \mathrm{d}\eta. \end{matrix} }[/math]

The ambient incident wavefield [math]\displaystyle{ \phi^{\mathrm{In}} }[/math] can also be expanded in the eigenfunctions corresponding to the incident wavefield upon [math]\displaystyle{ \Delta_l }[/math]. A detailed illustration of how to accomplish this will be given later. Let [math]\displaystyle{ D_{l0\mu}^{\mathrm{In}} }[/math] denote the coefficients of this ambient incident wavefield corresponding to the propagating modes and [math]\displaystyle{ D_{l\mu}^{\mathrm{In}} (\cdot) }[/math] denote the coefficients functions corresponding to the decaying modes (which are identically zero) of the incoming eigenfunction expansion for [math]\displaystyle{ \Delta_l }[/math]. The total incident wavefield upon body [math]\displaystyle{ \Delta_j }[/math] can now be expressed as

[math]\displaystyle{ \begin{matrix} &\phi_l^{\mathrm{I}}(r_l,\theta_l,z) = \phi^{\mathrm{In}}(r_l,\theta_l,z) + \sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \, \phi_j^{\mathrm{S}} (r_l,\theta_l,z)\\ &= \mathrm{e}^{\alpha z} \sum_{\mu = -\infty}^{\infty} \Big[ D_{l0\mu}^{\mathrm{In}} + \sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j H_{\nu-\mu}^{(1)} (\alpha R_{jl}) \mathrm{e}^{\i (\nu - \mu) \vartheta_{jl}} \Big] J_\mu (\alpha r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l}\\ & + \int\limits_0^{\infty} \psi (z,\eta) \sum_{\mu = -\infty}^{\infty} \Big[ D_{l\mu}^{\mathrm{In}}(\eta) + \sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \sum_{\nu = -\infty}^{\infty} A_{\nu}^j (\eta) (-1)^\mu K_{\nu - \mu} (\eta R_{jl}) \mathrm{e}^{\mathrm{i}(\nu - \mu) \vartheta_{jl}} \Big] I_\mu (\eta r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l} \mathrm{d}\eta. \end{matrix} }[/math]

The coefficients of the total incident potential upon [math]\displaystyle{ \Delta_l }[/math] are therefore given by \begin{subequations} (inc_coeff)

[math]\displaystyle{ \begin{matrix} D_{0\mu}^l &= D_{l0\mu}^{\mathrm{In}} + \sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j H_{\nu-\mu}^{(1)} (\alpha R_{jl}) \mathrm{e}^{\i (\nu - \mu) \vartheta_{jl}},\\ D_{\mu}^l(\eta) &= D_{l\mu}^{\mathrm{In}}(\eta) + \sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \sum_{\nu = -\infty}^{\infty} A_{\nu}^j (\eta) (-1)^\mu K_{\nu - \mu} (\eta R_{jl}) \mathrm{e}^{\mathrm{i}(\nu - \mu) \vartheta_{jl}}. \end{matrix} }[/math]

\end{subequations}

In general, it is possible to relate the total incident and scattered partial waves for any body through the diffraction characteristics of that body in isolation. There exist diffraction transfer operators [math]\displaystyle{ B_l }[/math] that relate the coefficients of the incident and scattered partial waves, such that

[math]\displaystyle{ (eq_B) A_l = B_l (D_l), \quad l=1, \ldots, N, }[/math]

where [math]\displaystyle{ A_l }[/math] are the scattered modes due to the incident modes [math]\displaystyle{ D_l }[/math]. In the case of a countable number of modes, (i.e. when the depth is finite), [math]\displaystyle{ B_l }[/math] is an infinite dimensional matrix. When the modes are functions of a continuous variable (i.e. infinite depth), [math]\displaystyle{ B_l }[/math] is the kernel of an integral operator. For the propagating and the decaying modes respectively, the scattered potential can be related by diffraction transfer operators acting in the following ways, \begin{subequations} (diff_op)

[math]\displaystyle{ \begin{matrix} A_{0\nu}^l &= \sum_{\mu = -\infty}^{\infty} B_{l\nu\mu}^\mathrm{pp} D_{0\mu}^l + \int\limits_{0}^{\infty} \sum_{\mu = -\infty}^{\infty} B_{l\nu\mu}^\mathrm{pd} (\xi) D_{\mu}^l (\xi) \mathrm{d}\xi,\\ A_\nu^l (\eta) &= \sum_{\mu = -\infty}^{\infty} B_{l\nu\mu}^\mathrm{dp} (\eta) D_{0\mu}^l + \int\limits_{0}^{\infty} \sum_{\mu = -\infty}^{\infty} B_{l\nu\mu}^\mathrm{dd} (\eta;\xi) D_{\mu}^l (\xi) \mathrm{d}\xi. \end{matrix} }[/math]

\end{subequations} The superscripts [math]\displaystyle{ \mathrm{p} }[/math] and [math]\displaystyle{ \mathrm{d} }[/math] are used to distinguish between propagating and decaying modes, the first superscript denotes the kind of scattered mode, the second one the kind of incident mode. If the diffraction transfer operators are known (their calculation will be discussed later), the substitution of equations (inc_coeff) into equations (diff_op) give the required equations to determine the coefficients and coefficient functions of the scattered wavefields of all bodies, \begin{subequations} (eq_op)

[math]\displaystyle{ \begin{matrix} &\begin{aligned} &A_{0n}^l = \sum_{\mu = -\infty}^{\infty} B_{ln\mu}^\mathrm{pp} \Big[ D_{l0\mu}^{\mathrm{In}} + \sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j H_{\nu-\mu}^{(1)} (\alpha R_{jl}) \mathrm{e}^{\i (\nu - \mu) \vartheta_{jl}} \Big]\\ & \ + \int\limits_{0}^{\infty} \sum_{\mu = -\infty}^{\infty} B_{ln\mu}^\mathrm{pd} (\xi) \Big[D_{l\mu}^{\mathrm{In}}(\eta) + \sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \sum_{\nu = -\infty}^{\infty} A_{\nu}^j (\eta) (-1)^\mu K_{\nu - \mu} (\eta R_{jl}) \mathrm{e}^{\mathrm{i}(\nu - \mu) \vartheta_{jl}} \Big] \mathrm{d}\xi, \end{aligned}\\ &\begin{aligned} &A_n^l (\eta) = \sum_{\mu = -\infty}^{\infty} B_{ln\mu}^\mathrm{dp} (\eta) \Big[ D_{l0\mu}^{\mathrm{In}} + \sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \sum_{\nu = - \infty}^{\infty} A_{0\nu}^j H_{\nu-\mu}^{(1)} (\alpha R_{jl}) \mathrm{e}^{\i (\nu - \mu) \vartheta_{jl}}\Big]\\ & \ + \int\limits_{0}^{\infty} \sum_{\mu = -\infty}^{\infty} B_{ln\mu}^\mathrm{dd} (\eta;\xi) \Big[ D_{l\mu}^{\mathrm{In}}(\eta) + \sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \sum_{\nu = -\infty}^{\infty} A_{\nu}^j (\eta) (-1)^\mu K_{\nu - \mu} (\eta R_{jl}) \mathrm{e}^{\mathrm{i}(\nu - \mu) \vartheta_{jl}}\Big] \mathrm{d}\xi, \end{aligned} \end{matrix} }[/math]

\end{subequations} [math]\displaystyle{ n \in \mathds{Z},\, l = 1, \ldots, N }[/math]. It has to be noted that all equations are coupled so that it is necessary to solve for all scattered coefficients and coefficient functions simultaneously.

For numerical calculations, the infinite sums have to be truncated and the integrals must be discretised. Implying a suitable truncation, the four different diffraction transfer operators can be represented by matrices which can be assembled in a big matrix [math]\displaystyle{ \mathbf{B}_l }[/math],

[math]\displaystyle{ \mathbf{B}_l = \left[ \begin{matrix}{cc} \mathbf{B}_l^{\mathrm{pp}} & \mathbf{B}_l^{\mathrm{pd}}\\ \mathbf{B}_l^{\mathrm{dp}} & \mathbf{B}_l^{\mathrm{dd}} \end{matrix} \right], }[/math]

the infinite depth diffraction transfer matrix. Truncating the coefficients accordingly, defining [math]\displaystyle{ {\bf a}^l }[/math] to be the vector of the coefficients of the scattered potential of body [math]\displaystyle{ \Delta_l }[/math], [math]\displaystyle{ \mathbf{d}_l^{\mathrm{In}} }[/math] to be the vector of coefficients of the ambient wavefield, and making use of a coordinate transformation matrix [math]\displaystyle{ {\bf T}_{jl} }[/math] given by \begin{subequations} (T_elem_deep)

[math]\displaystyle{ \begin{matrix} ({\bf T}_{jl})_{pq} &= H_{p-q}^{(1)}(\alpha R_{jl}) \, \mathrm{e}^{\mathrm{i}(p-q) \vartheta_{jl}}\\ =for the propagating modes, and= ({\bf T}_{jl})_{pq} &= (-1)^{q} \, K_{p-q} (\eta R_{jl}) \, \mathrm{e}^{\i (p-q) \vartheta_{jl}} \end{matrix} }[/math]

\end{subequations} for the decaying modes, a linear system of equations for the unknown coefficients follows from equations (eq_op),

[math]\displaystyle{ (eq_B_inf) {\bf a}_l = {\bf \hat{B}}_l \Big( {\bf d}_l^{\mathrm{In}} + \sum_{\genfrac{}{}{0pt}{}{j=1}{j \neq l}}^{N} \trans {\bf T}_{jl} \, {\bf a}_j \Big), \quad l=1, \ldots, N, }[/math]

where the left superscript [math]\displaystyle{ \mathrm{t} }[/math] indicates transposition. The matrix [math]\displaystyle{ {\bf \hat{B}}_l }[/math] denotes the infinite depth diffraction transfer matrix [math]\displaystyle{ {\bf B}_l }[/math] in which the elements associated with decaying scattered modes have been multiplied with the appropriate integration weights depending on the discretisation of the continuous variable.


\subsection{Calculation of the diffraction transfer matrix for bodies of arbitrary geometry}

Before we can apply the interaction theory we require the diffraction transfer matrices [math]\displaystyle{ \mathbf{B}_j }[/math] which relate the incident and the scattered potential for a body [math]\displaystyle{ \Delta_j }[/math] in isolation. The elements of the diffraction transfer matrix, [math]\displaystyle{ ({\bf B}_j)_{pq} }[/math], are the coefficients of the [math]\displaystyle{ p }[/math]th partial wave of the scattered potential due to a single unit-amplitude incident wave of mode [math]\displaystyle{ q }[/math] upon [math]\displaystyle{ \Delta_j }[/math].

While \citeauthor{kagemoto86}'s interaction theory was valid for bodies of arbitrary shape, they did not explain how to actually obtain the diffraction transfer matrices for bodies which did not have an axisymmetric geometry. This step was performed by goo90 who came up with an explicit method to calculate the diffraction transfer matrices for bodies of arbitrary geometry in the case of finite depth. Utilising a Green's function they used the standard method of transforming the single diffraction boundary-value problem to an integral equation for the source strength distribution function over the immersed surface of the body. However, the representation of the scattered potential which is obtained using this method is not automatically given in the cylindrical eigenfunction expansion. To obtain such cylindrical eigenfunction expansions of the potential goo90 used the representation of the free surface finite depth Green's function given by black75 and fenton78. \citeauthor{black75} and \citeauthor{fenton78}'s representation of the Green's function was based on applying Graf's addition theorem to the eigenfunction representation of the free surface finite depth Green's function given by john2. Their representation allowed the scattered potential to be represented in the eigenfunction expansion with the cylindrical coordinate system fixed at the point of the water surface above the mean centre position of the body.

It should be noted that, instead of using the source strength distribution function, it is also possible to consider an integral equation for the total potential and calculate the elements of the diffraction transfer matrix from the solution of this integral equation. An outline of this method for water of finite depth is given by kashiwagi00. We will present here a derivation of the diffraction transfer matrices for the case infinite depth based on a solution for the source strength distribution function. However, an equivalent derivation would be possible based on the solution for the total velocity potential.

To calculate the diffraction transfer matrix in infinite depth, we require the representation of the infinite depth free surface Green's function in cylindrical eigenfunctions,

[math]\displaystyle{ (green_inf) G(r,\theta,z;s,\varphi,c) &= \frac{\mathrm{i}\alpha}{2} \, \mathrm{e}^{\alpha (z+c)} \sum_{\nu=-\infty}^{\infty} H_\nu^{(1)}(\alpha r) J_\nu(\alpha s) \mathrm{e}^{\mathrm{i}\nu (\theta - \varphi)}\\ &\quad + \frac{1}{\pi^2} \int\limits_0^{\infty} \psi(z,\eta) \frac{\eta^2}{\eta^2+\alpha^2} \psi(c,\eta) \sum_{\nu=-\infty}^{\infty} K_\nu(\eta r) I_\nu(\eta s) \mathrm{e}^{\mathrm{i}\nu (\theta - \varphi)} \mathrm{d}\eta, }[/math]

[math]\displaystyle{ r \gt s }[/math], given by malte03.

We assume that we have represented the scattered potential in terms of the source strength distribution [math]\displaystyle{ \varsigma^j }[/math] so that the scattered potential can be written as

[math]\displaystyle{ (int_eq_1) \phi_j^\mathrm{S}(\mathbf{y}) = \int\limits_{\Gamma_j} G (\mathbf{y},\mathbf{\zeta}) \, \varsigma^j (\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta}, \quad \mathbf{y} \in D, }[/math]

where [math]\displaystyle{ D }[/math] is the volume occupied by the water and [math]\displaystyle{ \Gamma_j }[/math] is the immersed surface of body [math]\displaystyle{ \Delta_j }[/math]. The source strength distribution function [math]\displaystyle{ \varsigma^j }[/math] can be found by solving an integral equation. The integral equation is described in Weh_Lait and numerical methods for its solution are outlined in Sarp_Isa. Substituting the eigenfunction expansion of the Green's function

(green_inf) into  (int_eq_1), the scattered potential can

be written as

[math]\displaystyle{ \begin{matrix} &\phi_j^\mathrm{S}(r_j,\theta_j,z) = \mathrm{e}^{\alpha z} \sum_{\nu = - \infty}^{\infty} \bigg[ \frac{\mathrm{i}\alpha}{2} \int\limits_{\Gamma_j} \mathrm{e}^{\alpha c} J_\nu(\alpha s) \mathrm{e}^{-\mathrm{i}\nu \varphi} \varsigma^j(\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta} \bigg] H_\nu^{(1)} (\alpha r_j) \mathrm{e}^{\mathrm{i}\nu \theta_j}\\ & \, + \int\limits_{0}^{\infty} \psi(z,\eta) \sum_{\nu = - \infty}^{\infty} \bigg[ \frac{1}{\pi^2} \frac{\eta^2 }{\eta^2 + \alpha^2} \int\limits_{\Gamma_j} \psi(c,\eta) I_\nu(\eta s) \mathrm{e}^{-\mathrm{i}\nu \varphi} \varsigma^j({\bf{\zeta}}) \mathrm{d}\sigma_{\bf{\zeta}} \bigg] K_\nu (\eta r_j) \mathrm{e}^{\mathrm{i}\nu \theta_j} \mathrm{d}\eta, \end{matrix} }[/math]

where [math]\displaystyle{ \mathbf{\zeta}=(s,\varphi,c) }[/math] and [math]\displaystyle{ r\gt s }[/math]. This restriction implies that the eigenfunction expansion is only valid outside the escribed cylinder of the body.

The columns of the diffraction transfer matrix are the coefficients of the eigenfunction expansion of the scattered wavefield due to the different incident modes of unit-amplitude. The elements of the diffraction transfer matrix of a body of arbitrary shape are therefore given by \begin{subequations} (B_elem)

[math]\displaystyle{ \begin{matrix} ({\bf B}_j)_{pq} &= \frac{\mathrm{i}\alpha}{2} \int\limits_{\Gamma_j} \mathrm{e}^{\alpha c} J_p(\alpha s) \mathrm{e}^{-\mathrm{i}p \varphi} \varsigma_q^j(\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta}\\ =and= ({\bf B}_j)_{pq} &= \frac{1}{\pi^2} \frac{\eta^2}{\eta^2 + \alpha^2} \int\limits_{\Gamma_j} \psi(c,\eta) I_p(\eta s) \mathrm{e}^{-\mathrm{i}p \varphi} \varsigma_q^j(\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta} \end{matrix} }[/math]

\end{subequations} for the propagating and the decaying modes respectively, where [math]\displaystyle{ \varsigma_q^j(\mathbf{\zeta}) }[/math] is the source strength distribution due to an incident potential of mode [math]\displaystyle{ q }[/math] of the form \begin{subequations} (test_modes_inf)

[math]\displaystyle{ \begin{matrix} \phi_q^{\mathrm{I}}(s,\varphi,c) &= \mathrm{e}^{\alpha c} H_q^{(1)} (\alpha s) \mathrm{e}^{\mathrm{i}q \varphi}\\ =for the propagating modes, and= \phi_q^{\mathrm{I}}(s,\varphi,c) &= \psi(c,\eta) K_q (\eta s) \mathrm{e}^{\mathrm{i}q \varphi} \end{matrix} }[/math]

\end{subequations} for the decaying modes.

The diffraction transfer matrix of rotated bodies

For a non-axisymmetric body, a rotation about the mean centre position in the [math]\displaystyle{ (x,y) }[/math]-plane will result in a different diffraction transfer matrix. We will show how the diffraction transfer matrix of a body rotated by an angle [math]\displaystyle{ \beta }[/math] can be easily calculated from the diffraction transfer matrix of the non-rotated body. The rotation of the body influences the form of the elements of the diffraction transfer matrices in two ways. Firstly, the angular dependence in the integral over the immersed surface of the body is altered and, secondly, the source strength distribution function is different if the body is rotated. However, the source strength distribution function of the rotated body can be obtained by calculating the response of the non-rotated body due to rotated incident potentials. It will be shown that the additional angular dependence can be easily factored out of the elements of the diffraction transfer matrix.

The additional angular dependence caused by the rotation of the incident potential can be factored out of the normal derivative of the incident potential such that

[math]\displaystyle{ \frac{\partial \phi_{q\beta}^{\mathrm{I}}}{\partial n} = \frac{\partial \phi_{q}^{\mathrm{I}}}{\partial n} \mathrm{e}^{\mathrm{i}q \beta}, }[/math]

where [math]\displaystyle{ \phi_{q\beta}^{\mathrm{I}} }[/math] is the rotated incident potential. Since the integral equation for the determination of the source strength distribution function is linear, the source strength distribution function due to the rotated incident potential is thus just given by

[math]\displaystyle{ \varsigma_{q\beta}^j = \varsigma_q^j \, \mathrm{e}^{\mathrm{i}q \beta}. }[/math]

This is also the source strength distribution function of the rotated body due to the standard incident modes.

The elements of the diffraction transfer matrix [math]\displaystyle{ \mathbf{B}_j }[/math] are given by equations (B_elem). Keeping in mind that the body is rotated by the angle [math]\displaystyle{ \beta }[/math], the elements of the diffraction transfer matrix of the rotated body are given by \begin{subequations} (B_elem_rot)

[math]\displaystyle{ \begin{matrix} (\mathbf{B}_j^\beta)_{pq} &= \frac{\mathrm{i}\alpha}{2} \int\limits_{\Gamma_j} \mathrm{e}^{\alpha c} J_p(\alpha s) \mathrm{e}^{-\mathrm{i}p (\varphi+\beta)} \varsigma_{q\beta}^j(\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta},\\ =and= (\mathbf{B}_j^\beta)_{pq} &= \frac{1}{\pi^2} \frac{\eta^2}{\eta^2 + \alpha^2} \int\limits_{\Gamma_j} \psi(c,\eta) I_p(\eta s) \mathrm{e}^{-\mathrm{i}p (\varphi+\beta)} \varsigma_{q\beta}^j(\mathbf{\zeta}) \mathrm{d}\sigma_\mathbf{\zeta}, \end{matrix} }[/math]

\end{subequations} for the propagating and decaying modes respectively.

Thus the additional angular dependence caused by the rotation of the body can be factored out of the elements of the diffraction transfer matrix. The elements of the diffraction transfer matrix corresponding to the body rotated by the angle [math]\displaystyle{ \beta }[/math], [math]\displaystyle{ \mathbf{B}_j^\beta }[/math], are given by

[math]\displaystyle{ (B_rot) (\mathbf{B}_j^\beta)_{pq} = (\mathbf{B}_j)_{pq} \, \mathrm{e}^{\mathrm{i}(q-p) \beta}. }[/math]

As before, [math]\displaystyle{ (\mathbf{B})_{pq} }[/math] is understood to be the element of [math]\displaystyle{ \mathbf{B} }[/math] which corresponds to the coefficient of the [math]\displaystyle{ p }[/math]th scattered mode due to a unit-amplitude incident wave of mode [math]\displaystyle{ q }[/math]. Equation

(B_rot) applies to propagating and decaying modes likewise.


\subsection{Representation of the ambient wavefield in the eigenfunction representation} In Cartesian coordinates centred at the origin, the ambient wavefield is given by

[math]\displaystyle{ \phi^{\mathrm{In}}(x,y,z) = A \frac{g}{\omega} \, \mathrm{e}^{\mathrm{i}\alpha (x \cos \chi + y \sin \chi)+ \alpha z}, }[/math]

where [math]\displaystyle{ A }[/math] is the amplitude (in displacement) and [math]\displaystyle{ \chi }[/math] is the angle between the [math]\displaystyle{ x }[/math]-axis and the direction in which the wavefield travels. The interaction theory requires that the ambient wavefield, which is incident upon all bodies, is represented in the eigenfunction expansion of an incoming wave in the local coordinates of the body. The ambient wave can be represented in an eigenfunction expansion centred at the origin as

[math]\displaystyle{ \phi^{\mathrm{In}}(x,y,z) = A \frac{g}{\omega} \mathrm{e}^{\alpha z} \sum_{\mu = -\infty}^{\infty} \mathrm{e}^{\mathrm{i}\mu (\pi/2 - \theta + \chi)} J_\mu(\alpha r) }[/math]

\cite[p. 169]{linton01}. Since the local coordinates of the bodies are centred at their mean centre positions [math]\displaystyle{ O_l = (O_x^l,O_y^l) }[/math], a phase factor has to be defined which accounts for the position from the origin. Including this phase factor the ambient wavefield at the [math]\displaystyle{ l }[/math]th body is given by

[math]\displaystyle{ \phi^{\mathrm{In}}(r_l,\theta_l,z) = A \frac{g}{\omega} \, e^{\mathrm{i}\alpha (O_x^l \cos \chi + O_x^l \sin \chi)} \, \mathrm{e}^{\alpha z} \sum_{\mu = -\infty}^{\infty} \mathrm{e}^{\mathrm{i}\mu (\pi/2 - \chi)} J_\mu(\alpha r_l) \mathrm{e}^{\mathrm{i}\mu \theta_l}. }[/math]

Solving the resulting system of equations

After the coefficient vector of the ambient incident wavefield, the diffraction transfer matrices and the coordinate transformation matrices have been calculated, the system of equations (eq_B_inf), has to be solved. This system can be represented by the following matrix equation,

[math]\displaystyle{ \left[ \begin{matrix}{c} {\bf a}_1\\ {\bf a}_2\\ \\ \vdots \\ \\ {\bf a}_N \end{matrix} \right] = \left[ \begin{matrix}{c} \hat{{\bf B}}_1 {\bf d}_1^\mathrm{In}\\ \hat{{\bf B}}_2 {\bf d}_2^\mathrm{In}\\ \\ \vdots \\ \\ \hat{{\bf B}}_N {\bf d}_N^\mathrm{In} \end{matrix} \right]+ \left[ \begin{matrix}{ccccc} \mathbf{0} & \hat{{\bf B}}_1 \trans {\bf T}_{21} & \hat{{\bf B}}_1 \trans {\bf T}_{31} & \dots & \hat{{\bf B}}_1 \trans {\bf T}_{N1}\\ \hat{{\bf B}}_2 \trans {\bf T}_{12} & \mathbf{0} & \hat{{\bf B}}_2 \trans {\bf T}_{32} & \dots & \hat{{\bf B}}_2 \trans {\bf T}_{N2}\\ & & \mathbf{0} & &\\ \vdots & & & \ddots & \vdots\\ & & & & \\ \hat{{\bf B}}_N \trans {\bf T}_{1N} & & \dots & & \mathbf{0} \end{matrix} \right] \left[ \begin{matrix}{c} {\bf a}_1\\ {\bf a}_2\\ \\ \vdots \\ \\ {\bf a}_N \end{matrix} \right], }[/math]

where [math]\displaystyle{ \mathbf{0} }[/math] denotes the zero-matrix which is of the same dimension as [math]\displaystyle{ \hat{{\bf B}}_j }[/math], say [math]\displaystyle{ n }[/math]. This matrix equation can be easily transformed into a classical [math]\displaystyle{ (N \, n) }[/math]-dimensional linear system of equations.