Kagemoto and Yue Interaction Theory: Difference between revisions

From WikiWaves
Jump to navigationJump to search
No edit summary
 
(122 intermediate revisions by 5 users not shown)
Line 1: Line 1:
= Introduction =
This is an interaction theory which provides the exact solution (i.e. it is not based on a [[Wide Spacing Approximation]]).
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 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]]. [[Interaction Theory for Cylinders]] presents a simplified version for cylinders.  


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 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 coordinates 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
The theory is described in [[Kagemoto and Yue 1986]] and in
[[Peter and Meylan 2004]].
[[Peter and Meylan 2004]].
 
The derivation of the theory in [[Infinite Depth]] is also presented, see
[[Kagemoto and Yue Interaction Theory for Infinite Depth]].
   
   
[[Category:Linear Water-Wave Theory]]
[[Category:Interaction Theory]]
 
 


= Equations of Motion =


The problem consists of <math>n</math> bodies
<math>\Delta_j</math> with immersed body
surface <math>\Gamma_j</math>. Each body is subject to
the [[Standard Linear Wave Scattering Problem]] and the particluar
equations of motion for each body (e.g. rigid, or freely floating)
can be different for each body.
It is a [[Frequency Domain Problem]] with frequency <math>\omega</math>.
The solution is exact, up to the
restriction that the escribed cylinder of each body may not contain any
other body.
To simplify notation, <math>\mathbf{y} = (x,y,z)</math> always denotes a point
in the water, which is assumed to be of [[Finite Depth]] <math>h</math>,
while <math>\mathbf{x}</math> always denotes a point of the undisturbed water
surface assumed at <math>z=0</math>.


{{standard linear wave scattering equations}}


We extend the finite depth interaction theory of [[kagemoto86]] to
The [[Sommerfeld Radiation Condition]] is also imposed.
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.


=Eigenfunction expansion of the potential=


 
Each body is subject to an incident potential and moves in response to this
==Introduction==
incident potential to produce a scattered potential. Each of these is
The scattering of water waves by floating or submerged
expanded using the [[Cylindrical Eigenfunction Expansion]]
bodies is of wide practical importance.
The scattered potential of a body
Although the problem is non-linear, if the
<math>\Delta_j</math> can be expressed as
wave amplitude is sufficiently small, the
<center><math>  
problem can be linearised. The linear problem is still the basis of
\phi_j^\mathrm{S} (r_j,\theta_j,z) =  
the engineering design of most off-shore structures and is the standard
\sum_{m=0}^{\infty} f_m(z) \sum_{\mu = -
model of geophysical phenomena such as the wave forcing of ice floes.
\infty}^{\infty} A_{m \mu}^j K_\mu (k_m r_j) \mathrm{e}^{\mathrm{i}\mu \theta_j},
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.
 
 
 
==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>
</math></center>
with discrete coefficients <math>A_{m \nu}^j</math>. The positive wavenumber <math>k</math>
with discrete coefficients <math>A_{m \mu}^j</math>, where <math>(r_j,\theta_j,z)</math>
is related to <math>\alpha</math> by the dispersion relation
are cylindrical polar coordinates centered at each body
<center><math> (eq_k)
<center><math>
\alpha = k \tanh k d,
f_m(z) = \frac{\cos k_m (z+h)}{\cos k_m h}.
</math></center>
</math></center>
and the values of <math>k_m</math>, <math>m>0</math>, are given as positive real roots of
where <math>k_m</math> are found from <math>\alpha</math> by the [[Dispersion Relation for a Free Surface]]
the dispersion relation
<center><math>  
<center><math> (eq_k_m)
\alpha + k_m \tan k_m h = 0\,.
\alpha + k_m \tan k_m d = 0.
</math></center>  
</math></center>  
where <math>k_0</math> is the
imaginary root with negative imaginary part
and <math>k_m</math>, <math>m>0</math>, are given the positive real roots ordered
with increasing size.
The incident potential upon body <math>\Delta_j</math> can be also be expanded in
The incident potential upon body <math>\Delta_j</math> can be also be expanded in
cylindrical eigenfunctions,  
regular cylindrical eigenfunctions,  
<center><math> (basisrep_in_d)
<center><math>  
\phi_j^\mathrm{I} (r_j,\theta_j,z) &= \frac{\cosh k(z+d)}{\cosh kd}
\phi_j^\mathrm{I} (r_j,\theta_j,z) = \sum_{n=0}^{\infty} f_n(z)
\sum_{\mu = - \infty}^{\infty} D_{0 \mu}^j J_\mu (k r_j) \mathrm{e}^{\mathrm{i}\mu
\sum_{\nu = - \infty}^{\infty} D_{n\nu}^j I_\nu (k_n r_j) \mathrm{e}^{\mathrm{i}\nu \theta_j},
\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>
</math></center>
given by [[black75]] and [[fenton78]], needs to be used instead
with discrete coefficients <math>D_{n\nu}^j</math>. In these expansions, <math>I_\nu</math>
of the infinite depth Green's function  (green_inf).
and <math>K_\nu</math> denote the modified [http://en.wikipedia.org/wiki/Bessel_function : Bessel functions]
The elements of <math>{\bf B}_j</math> are therefore given by
of the first and second kind, respectively, both of order <math>\nu</math>.
\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.
 
 


Note that the term for <math>m =0</math> or
<math>n=0</math> corresponds to the propagating modes while the
terms for <math>m\geq 1</math> (<math>n\geq 1</math>) correspond to the evanescent modes.


==Wave forcing of an ice floe of arbitrary geometry==
=Derivation of the system of equations=


The interaction theory which has been developed so far has been
A system of equations for the unknown
for arbitrary bodies. No assumption has been made about the body
coefficients of the
geometry or its equations of motion. However, we will now use this
scattered wavefields of all bodies is developed. This system of
interaction theory to make calculations for the specific case of ice
equations is based on transforming the
floes. Ice floes form in vast fields consisting of hundreds if not
scattered potential of <math>\Delta_j</math> into an incident potential upon
thousands of individual floes and furthermore  most ice floe fields
<math>\Delta_l</math> (<math>j \neq l</math>). Doing this for all bodies simultaneously,
occur in the deep ocean. For this reason they are ideally suited to
and relating the incident and scattered potential for each body, a system
the application of the scattering theory we have just developed.  
of equations for the unknown coefficients is developed.  
Furthermore, the presence of the ice lengthens the wavelength
Making use of the periodicity of the geometry and of the ambient incident
making it more difficult to determine how deep the water must
wave, this system of equations can then be simplified.
be to be approximately infinite.


===Mathematical model for ice floes===
The scattered potential <math>\phi_j^{\mathrm{S}}</math> of body <math>\Delta_j</math> needs to be
We will briefly describe the mathematical model which is used to
represented in terms of the incident potential <math>\phi_l^{\mathrm{I}}</math>
describe ice floes. A more detailed account can be found in
upon <math>\Delta_l</math>, <math>j \neq l</math>. This can be accomplished by using
[[Squire_review]]. We assume that the ice floe is sufficiently thin
[[Graf's Addition Theorem]]
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>
<center><math>
\frac{\partial^2 w}{\partial n^2} + \nu \frac{\partial^2 w}{\partial
K_\tau(k_m r_j) \mathrm{e}^{\mathrm{i}\tau (\theta_j-\varphi_{jl})} =
s^2} = 0 \quad =and=  \quad \frac{\partial^3 w}{\partial n^3} + (2 - \nu)
\sum_{\nu = - \infty}^{\infty} K_{\tau + \nu} (k_m R_{jl}) \,
\frac{\partial^3 w}{\partial n \partial s^2} = 0, \quad
I_\nu (k_m r_l) \mathrm{e}^{\mathrm{i}\nu (\pi - \theta_l + \varphi_{jl})}, \quad j \neq l,
\mathbf{x} \in \partial \Delta,
</math></center>
</math></center>
where <math>n</math> and <math>s</math> denote the normal and tangential directions on
which is valid provided that <math>r_l < R_{jl}</math>. Here, <math>(R_{jl},\varphi_{jl})</math> are the polar coordinates of the mean centre position of <math>\Delta_{l}</math> in the local coordinates of <math>\Delta_{j}</math>.  
<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,
The limitation <math>r_l < R_{jl}</math> only requires that the escribed cylinder of each body
<center><math>
<math>\Delta_l</math> does not enclose any other origin <math>O_j</math> (<math>j \neq l</math>). However, the
(\bar{x},\bar{y},\bar{z}) = \frac{1}{a} (x,y,z), \quad \bar{w} =
expansion of the scattered and incident potential in cylindrical
\frac{w}{a}, \quad \bar{\alpha} = a\, \alpha, \quad \bar{\omega} = \omega
eigenfunctions is only valid outside the escribed cylinder of each
\sqrt{\frac{a}{g}} \quad =and=  \quad \bar{\phi} = \frac{\phi}{a
body. Therefore the condition that the
\sqrt{a g}},
escribed cylinder of each body <math>\Delta_l</math> does not enclose any other
</math></center>
origin <math>O_j</math> (<math>j \neq l</math>) is superseded by the more rigorous
where <math>a</math> is a length parameter associated with the floe.
restriction that the escribed cylinder of each body may not contain any
In non-dimensional variables, the equation for the ice floe
other body.  
(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.
Making use of the eigenfunction expansion as well as [[Graf's Addition Theorem]], the scattered potential
The water velocity potential must satisfy the boundary value problem
of <math>\Delta_j</math> can be expressed in terms of the
\begin{subequations} (water)
incident potential upon <math>\Delta_l</math> as
<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>
<center><math>
\lim_{\tilde{r} \rightarrow \infty} \sqrt{\tilde{r}} \, \Big(
\phi_j^{\mathrm{S}} (r_l,\theta_l,z)
\frac{\partial}{\partial \tilde{r}} - \mathrm{i}k
= \sum_{m=0}^\infty f_m(z) \sum_{\tau = -
\Big) (\phi - \phi^{\mathrm{In}}) = 0,
\infty}^{\infty} A_{m\tau}^j \sum_{\nu = -\infty}^{\infty}
(-1)^\nu K_{\tau-\nu} (k_m  R_{jl}) I_\nu (k_m r_l) \mathrm{e}^{\mathrm{i}\nu
\theta_l} \mathrm{e}^{\mathrm{i}(\tau-\nu) \varphi_{jl}}  
</math></center>
</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>
<center><math>
\frac{\partial \phi}{\partial z} = 0, \quad \mathbf{y} \in D,\: z=-d.
= \sum_{m=0}^\infty f_m(z) \sum_{\nu =
-\infty}^{\infty} \Big[ \sum_{\tau = - \infty}^{\infty} A_{m\tau}^j
(-1)^\nu K_{\tau-\nu} (k_m R_{jl}) \mathrm{e}^{\mathrm{i}(\tau - \nu)
\varphi_{jl}} \Big] I_\nu (k_m r_l) \mathrm{e}^{\mathrm{i}\nu \theta_l}.
</math></center>
</math></center>
In water of finite depth, the positive real wavenumber <math>k</math> is related
The ambient incident wavefield <math>\phi^{\mathrm{In}}</math> can also be
to the radian frequency by the dispersion relations  (eq_k).
expanded in the eigenfunctions corresponding to the incident wavefield upon
 
<math>\Delta_l</math>. Let <math>\tilde{D}_{n\nu}^{l}</math> denote the coefficients of this
 
ambient incident wavefield in the incoming eigenfunction expansion for
===The wavelength under the ice floe=== (sec:kappa)
<math>\Delta_l</math> (cf. the example in [[Cylindrical Eigenfunction Expansion]]).  
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>
<center><math>
\kappa^* \tan \kappa^* d = - \frac{\alpha}{\beta \kappa^{*4} - \gamma
\phi^{\mathrm{In}}(r_l,\theta_l,z)= \sum_{n=0}^\infty f_n(z) \sum_{\nu = -\infty}^{\infty}
\alpha +1},
\tilde{D}_{n\nu}^{l}  I_\nu (k_n
r_l) \mathrm{e}^{\mathrm{i}\nu \theta_l}.
</math></center>
</math></center>
where <math>\kappa^*</math> is the wavenumber under the plate. The purely imaginary
The total
roots of this dispersion relation correspond to the propagating modes
incident wavefield upon body <math>\Delta_j</math> can now be expressed as  
and their absolute value is given as the positive root of
<center><math>
<center><math>
\kappa \tanh \kappa d = \frac{\alpha}{\beta \kappa^4 - \gamma
\phi_l^{\mathrm{I}}(r_l,\theta_l,z) = \phi^{\mathrm{In}}(r_l,\theta_l,z) +
\alpha + 1}.
\sum_{j=1,j \neq l}^{N} \, \phi_j^{\mathrm{S}}
</math></center>
(r_l,\theta_l,z)
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>
</math></center>
In the case of a shallow draft, the fact that the Green's function is
This allows us to write
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>
<center><math>
\nabla^4 w^k = \lambda_k w^k.
\sum_{n=0}^{\infty} f_n(z)
\sum_{\nu = - \infty}^{\infty} D_{n\nu}^l I_\nu (k_n r_l) \mathrm{e}^{\mathrm{i}\nu \theta_l}
</math></center>
</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>
<center><math>
\int\limits_\Delta w^j (\mathbf{\xi}) w^k (\mathbf{\xi})
= \sum_{n=0}^\infty f_n(z) \sum_{\nu = -\infty}^{\infty}
\mathrm{d}\sigma_{\mathbf{\xi}} = \delta _{jk},
  \Big[  \tilde{D}_{n\nu}^{l} +
</math></center>
\sum_{j=1,j \neq  l}^{N} \sum_{\tau =
where <math>\delta _{jk}</math> is the Kronecker delta. The eigenvalues <math>\lambda_k</math>
-\infty}^{\infty} A_{n\tau}^j (-1)^\nu K_{\tau - \nu} (k_n
have the property that <math>\lambda_k \rightarrow \infty</math> as </math>k \rightarrow
R_{jl})  \mathrm{e}^{\mathrm{i}(\tau - \nu) \varphi_{jl}} \Big] I_\nu (k_n
\infty<math> and we order the modes by increasing eigenvalue. These modes can be
r_l) \mathrm{e}^{\mathrm{i}\nu \theta_l}.
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>
</math></center>
which is a matrix equation in <math>c_k</math>.
It therefore follows that
 
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>
<center><math>
E_2 = \left( \, \int\limits_{\Delta}
D_{n\nu}^l  =  
\big| w_{i}(\mathbf{x})-w_{f}(\mathbf{x}) \big|^2 \mathrm{d}\mathbf{x} \,
  \tilde{D}_{n\nu}^{l} +
\right)^{1/2},
\sum_{j=1,j \neq  l}^{N} \sum_{\tau =
-\infty}^{\infty} A_{n\tau}^j (-1)^\nu K_{\tau - \nu} (k_n
R_{jl}) \mathrm{e}^{\mathrm{i}(\tau - \nu) \varphi_{jl}}  
</math></center>
</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
= Final Equations =
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
The scattered and incident potential of each body <math>\Delta_l</math> can be related by the
now be considered for both finite and infinite depth.
[[Diffraction Transfer Matrix]] acting in the following way,
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>
<center><math>
\phi (r,\theta,z) &= \mathrm{e}^{\alpha z} \sum_{\nu = -
A_{m \mu}^l = \sum_{n=0}^{\infty} \sum_{\nu = -\infty}^{\infty} B_{m n
T_H}^{T_H} A_{0\nu} H_\nu^{(1)} (\alpha r) \mathrm{e}^{\mathrm{i}\nu
\mu \nu}^l D_{n\nu}^l.
\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>  
</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 =
The substitution of this into the equation for relating
\eta_2-\eta_1<math> as well as </math>h_{_{T_R}} =
the coefficients <math>D_{n\nu}^l</math> and
\eta_{_{T_R}}-\eta_{_{T_R-1}}</math>, which corresponds to the mid-point
<math>A_{m \mu}^l</math>gives the
quadrature rule.
required equations to determine the coefficients of the scattered
Different quadrature rules such as Gaussian quadrature
wavefields of all bodies,  
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>
<center><math>
\phi (r,\theta,z) &= \frac{\cosh k (z+d)}{\cosh kd} \sum_{\nu = -
A_{m\mu}^l = \sum_{n=0}^{\infty}
T_H}^{T_H} A_{0\nu} H_\nu^{(1)} (k r) \mathrm{e}^{\mathrm{i}\nu \theta}\\
\sum_{\nu = -\infty}^{\infty} B_{mn\mu\nu}^l
& \quad + \sum_{m = 1}^{T_R} \frac{\cos k_m(z+d)}{\cos k_m d}
\Big[ \tilde{D}_{n\nu}^{l} +
\sum_{\nu = - T_K}^{T_K} A_{m\nu} K_\nu (k_m r) \mathrm{e}^{\mathrm{i}\nu \theta}.
\sum_{j=1,j \neq  l}^{N} \sum_{\tau =
-\infty}^{\infty} A_{n\tau}^j (-1)^\nu K_{\tau - \nu(k_n
R_{jl}) \mathrm{e}^{\mathrm{i}(\tau -\nu) \varphi_{jl}} \Big],
</math></center>
</math></center>
In both cases, the dimension of the diffraction transfer matrix,
<math>m \in \mathbb{N}</math>, <math>\mu \in \mathbb{Z}</math>, <math>l=1,\dots,N</math>.
<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>

Latest revision as of 10:24, 2 May 2010

Introduction

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. Interaction Theory for Cylinders presents a simplified version for cylinders.

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

The derivation of the theory in Infinite Depth is also presented, see Kagemoto and Yue Interaction Theory for Infinite Depth.

Equations of Motion

The problem consists of [math]\displaystyle{ n }[/math] bodies [math]\displaystyle{ \Delta_j }[/math] with immersed body surface [math]\displaystyle{ \Gamma_j }[/math]. Each body is subject to the Standard Linear Wave Scattering Problem and the particluar equations of motion for each body (e.g. rigid, or freely floating) can be different for each body. It is a Frequency Domain Problem with frequency [math]\displaystyle{ \omega }[/math]. The solution is exact, up to the restriction that the escribed cylinder of each body may not contain any other body. To simplify notation, [math]\displaystyle{ \mathbf{y} = (x,y,z) }[/math] always denotes a point in the water, which is assumed to be of Finite Depth [math]\displaystyle{ h }[/math], while [math]\displaystyle{ \mathbf{x} }[/math] always denotes a point of the undisturbed water surface assumed at [math]\displaystyle{ z=0 }[/math].

The equations are the following

[math]\displaystyle{ \begin{align} \Delta\phi &=0, &-h\lt z\lt 0,\,\,\mathbf{x} \in \Omega \\ \partial_z\phi &= 0, &z=-h, \\ \partial_z \phi &= \alpha \phi, &z=0,\,\,\mathbf{x} \in \partial \Omega_{\mathrm{F}}, \end{align} }[/math]


(note that the last expression can be obtained from combining the expressions:

[math]\displaystyle{ \begin{align} \partial_z \phi &= -\mathrm{i} \omega \zeta, &z=0,\,\,\mathbf{x} \in \partial \Omega_{\mathrm{F}}, \\ \mathrm{i} \omega \phi &= g\zeta, &z=0,\,\,\mathbf{x} \in \partial \Omega_{\mathrm{F}}, \end{align} }[/math]

where [math]\displaystyle{ \alpha = \omega^2/g \, }[/math])

[math]\displaystyle{ \partial_n\phi = \mathcal{L}\phi, \quad \mathbf{x}\in\partial\Omega_B, }[/math]

where [math]\displaystyle{ \mathcal{L} }[/math] is a linear operator which relates the normal and potential on the body surface through the physics of the body.

The Sommerfeld Radiation Condition is also imposed.

Eigenfunction expansion of the potential

Each body is subject to an incident potential and moves in response to this incident potential to produce a scattered potential. Each of these is expanded using the Cylindrical Eigenfunction Expansion The scattered potential of a body [math]\displaystyle{ \Delta_j }[/math] can be expressed as

[math]\displaystyle{ \phi_j^\mathrm{S} (r_j,\theta_j,z) = \sum_{m=0}^{\infty} f_m(z) \sum_{\mu = - \infty}^{\infty} A_{m \mu}^j K_\mu (k_m r_j) \mathrm{e}^{\mathrm{i}\mu \theta_j}, }[/math]

with discrete coefficients [math]\displaystyle{ A_{m \mu}^j }[/math], where [math]\displaystyle{ (r_j,\theta_j,z) }[/math] are cylindrical polar coordinates centered at each body

[math]\displaystyle{ f_m(z) = \frac{\cos k_m (z+h)}{\cos k_m h}. }[/math]

where [math]\displaystyle{ k_m }[/math] are found from [math]\displaystyle{ \alpha }[/math] by the Dispersion Relation for a Free Surface

[math]\displaystyle{ \alpha + k_m \tan k_m h = 0\,. }[/math]

where [math]\displaystyle{ k_0 }[/math] is the imaginary root with negative imaginary part and [math]\displaystyle{ k_m }[/math], [math]\displaystyle{ m\gt 0 }[/math], are given the positive real roots ordered with increasing size.

The incident potential upon body [math]\displaystyle{ \Delta_j }[/math] can be also be expanded in regular cylindrical eigenfunctions,

[math]\displaystyle{ \phi_j^\mathrm{I} (r_j,\theta_j,z) = \sum_{n=0}^{\infty} f_n(z) \sum_{\nu = - \infty}^{\infty} D_{n\nu}^j I_\nu (k_n r_j) \mathrm{e}^{\mathrm{i}\nu \theta_j}, }[/math]

with discrete coefficients [math]\displaystyle{ D_{n\nu}^j }[/math]. In these expansions, [math]\displaystyle{ I_\nu }[/math] and [math]\displaystyle{ K_\nu }[/math] denote the modified : Bessel functions of the first and second kind, respectively, both of order [math]\displaystyle{ \nu }[/math].

Note that the term for [math]\displaystyle{ m =0 }[/math] or [math]\displaystyle{ n=0 }[/math] corresponds to the propagating modes while the terms for [math]\displaystyle{ m\geq 1 }[/math] ([math]\displaystyle{ n\geq 1 }[/math]) correspond to the evanescent modes.

Derivation of the system of equations

A system of equations for the unknown coefficients of the scattered wavefields of all bodies is 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 is developed. Making use of the periodicity of the geometry and of the ambient incident wave, this system of equations can then be simplified.

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]. This can be accomplished by using Graf's Addition Theorem

[math]\displaystyle{ K_\tau(k_m r_j) \mathrm{e}^{\mathrm{i}\tau (\theta_j-\varphi_{jl})} = \sum_{\nu = - \infty}^{\infty} K_{\tau + \nu} (k_m R_{jl}) \, I_\nu (k_m r_l) \mathrm{e}^{\mathrm{i}\nu (\pi - \theta_l + \varphi_{jl})}, \quad j \neq l, }[/math]

which is valid provided that [math]\displaystyle{ r_l \lt R_{jl} }[/math]. Here, [math]\displaystyle{ (R_{jl},\varphi_{jl}) }[/math] are the polar coordinates of the mean centre position of [math]\displaystyle{ \Delta_{l} }[/math] in the local coordinates of [math]\displaystyle{ \Delta_{j} }[/math].

The limitation [math]\displaystyle{ r_l \lt R_{jl} }[/math] 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 eigenfunction expansion as well as Graf's Addition Theorem, the scattered potential of [math]\displaystyle{ \Delta_j }[/math] can be expressed in terms of the incident potential upon [math]\displaystyle{ \Delta_l }[/math] as

[math]\displaystyle{ \phi_j^{\mathrm{S}} (r_l,\theta_l,z) = \sum_{m=0}^\infty f_m(z) \sum_{\tau = - \infty}^{\infty} A_{m\tau}^j \sum_{\nu = -\infty}^{\infty} (-1)^\nu K_{\tau-\nu} (k_m R_{jl}) I_\nu (k_m r_l) \mathrm{e}^{\mathrm{i}\nu \theta_l} \mathrm{e}^{\mathrm{i}(\tau-\nu) \varphi_{jl}} }[/math]
[math]\displaystyle{ = \sum_{m=0}^\infty f_m(z) \sum_{\nu = -\infty}^{\infty} \Big[ \sum_{\tau = - \infty}^{\infty} A_{m\tau}^j (-1)^\nu K_{\tau-\nu} (k_m R_{jl}) \mathrm{e}^{\mathrm{i}(\tau - \nu) \varphi_{jl}} \Big] I_\nu (k_m r_l) \mathrm{e}^{\mathrm{i}\nu \theta_l}. }[/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]. Let [math]\displaystyle{ \tilde{D}_{n\nu}^{l} }[/math] denote the coefficients of this ambient incident wavefield in the incoming eigenfunction expansion for [math]\displaystyle{ \Delta_l }[/math] (cf. the example in Cylindrical Eigenfunction Expansion).

[math]\displaystyle{ \phi^{\mathrm{In}}(r_l,\theta_l,z)= \sum_{n=0}^\infty f_n(z) \sum_{\nu = -\infty}^{\infty} \tilde{D}_{n\nu}^{l} I_\nu (k_n r_l) \mathrm{e}^{\mathrm{i}\nu \theta_l}. }[/math]

The total incident wavefield upon body [math]\displaystyle{ \Delta_j }[/math] can now be expressed as

[math]\displaystyle{ \phi_l^{\mathrm{I}}(r_l,\theta_l,z) = \phi^{\mathrm{In}}(r_l,\theta_l,z) + \sum_{j=1,j \neq l}^{N} \, \phi_j^{\mathrm{S}} (r_l,\theta_l,z) }[/math]

This allows us to write

[math]\displaystyle{ \sum_{n=0}^{\infty} f_n(z) \sum_{\nu = - \infty}^{\infty} D_{n\nu}^l I_\nu (k_n r_l) \mathrm{e}^{\mathrm{i}\nu \theta_l} }[/math]
[math]\displaystyle{ = \sum_{n=0}^\infty f_n(z) \sum_{\nu = -\infty}^{\infty} \Big[ \tilde{D}_{n\nu}^{l} + \sum_{j=1,j \neq l}^{N} \sum_{\tau = -\infty}^{\infty} A_{n\tau}^j (-1)^\nu K_{\tau - \nu} (k_n R_{jl}) \mathrm{e}^{\mathrm{i}(\tau - \nu) \varphi_{jl}} \Big] I_\nu (k_n r_l) \mathrm{e}^{\mathrm{i}\nu \theta_l}. }[/math]

It therefore follows that

[math]\displaystyle{ D_{n\nu}^l = \tilde{D}_{n\nu}^{l} + \sum_{j=1,j \neq l}^{N} \sum_{\tau = -\infty}^{\infty} A_{n\tau}^j (-1)^\nu K_{\tau - \nu} (k_n R_{jl}) \mathrm{e}^{\mathrm{i}(\tau - \nu) \varphi_{jl}} }[/math]

Final Equations

The scattered and incident potential of each body [math]\displaystyle{ \Delta_l }[/math] can be related by the Diffraction Transfer Matrix acting in the following way,

[math]\displaystyle{ A_{m \mu}^l = \sum_{n=0}^{\infty} \sum_{\nu = -\infty}^{\infty} B_{m n \mu \nu}^l D_{n\nu}^l. }[/math]

The substitution of this into the equation for relating the coefficients [math]\displaystyle{ D_{n\nu}^l }[/math] and [math]\displaystyle{ A_{m \mu}^l }[/math]gives the required equations to determine the coefficients of the scattered wavefields of all bodies,

[math]\displaystyle{ A_{m\mu}^l = \sum_{n=0}^{\infty} \sum_{\nu = -\infty}^{\infty} B_{mn\mu\nu}^l \Big[ \tilde{D}_{n\nu}^{l} + \sum_{j=1,j \neq l}^{N} \sum_{\tau = -\infty}^{\infty} A_{n\tau}^j (-1)^\nu K_{\tau - \nu} (k_n R_{jl}) \mathrm{e}^{\mathrm{i}(\tau -\nu) \varphi_{jl}} \Big], }[/math]

[math]\displaystyle{ m \in \mathbb{N} }[/math], [math]\displaystyle{ \mu \in \mathbb{Z} }[/math], [math]\displaystyle{ l=1,\dots,N }[/math].