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.
\
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)} \d\eta,
\lt math\gt 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})
\d\sigma_\mathbf{\zeta}, \quad \mathbf{y} \in D,
where \lt math\gt 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})
\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}})
\d\sigma_{\bf{\zeta}} \bigg] K_\nu (\eta r_j) \mathrm{e}^{\mathrm{i}\nu \theta_j} \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})
\d\sigma_\mathbf{\zeta}\\
\intertext{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}) \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}\\
\intertext{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},
where \lt math\gt \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}.
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 \lt math\gt \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}) \d\sigma_\mathbf{\zeta},\\
\intertext{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}) \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}.
As before, \lt math\gt (\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},
where \lt math\gt 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)
\cite[p. 169]{linton01}.
Since the local coordinates of the bodies are centred at their mean
centre positions \lt math\gt 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}.
===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,
\lt center\gt \lt math\gt
\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 \lt math\gt \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.
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]\displaystyle{ d }[/math], the scattered potential of a body
[math]\displaystyle{ \Delta_j }[/math] can be expanded in cylindrical eigenfunctions,
[math]\displaystyle{ (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 \lt math\gt A_{m \nu}^j }[/math]. The positive wavenumber [math]\displaystyle{ k }[/math]
is related to [math]\displaystyle{ \alpha }[/math] by the dispersion relation
[math]\displaystyle{ (eq_k)
\alpha = k \tanh k d,
and the values of \lt math\gt k_m }[/math], [math]\displaystyle{ m\gt 0 }[/math], are given as positive real roots of
the dispersion relation
[math]\displaystyle{ (eq_k_m)
\alpha + k_m \tan k_m d = 0.
The incident potential upon body \lt math\gt \Delta_j }[/math] can be also be expanded in
cylindrical eigenfunctions,
[math]\displaystyle{ (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 \lt math\gt 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
[math]\displaystyle{
{\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 \lt math\gt {\bf a}_l }[/math] is the coefficient vector of the scattered
wave, [math]\displaystyle{ {\bf d}_l^\mathrm{In} }[/math] is the coefficient vector of the
ambient incident wave, [math]\displaystyle{ {\bf B}_l }[/math] is the diffraction transfer
matrix of [math]\displaystyle{ \Delta_l }[/math] and [math]\displaystyle{ {\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
[math]\displaystyle{ (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 \lt math\gt {\bf B}_j }[/math] are therefore given by
\begin{subequations} (B_elem_d)
[math]\displaystyle{ \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})
\d\sigma_\mathbf{\zeta}\\
\intertext{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}) \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_d)
[math]\displaystyle{ \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}\\
\intertext{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]
\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]\displaystyle{ W }[/math] is that of the water surface and [math]\displaystyle{ W }[/math] is required to satisfy the linear
plate equation in the area occupied by the ice floe [math]\displaystyle{ \Delta }[/math]. In analogy to
(time), [math]\displaystyle{ 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
[math]\displaystyle{ (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 \lt math\gt \rho }[/math], the modulus of rigidity of the
ice floe [math]\displaystyle{ D }[/math], its density [math]\displaystyle{ \rho_\Delta }[/math] and its
thickness [math]\displaystyle{ h }[/math]. The right-hand-side of (plate_non) arises from the
linearised Bernoulli equation. It needs to be recalled that
[math]\displaystyle{ \mathbf{x} }[/math] always denotes a point of the undisturbed water surface.
Free edge boundary conditions apply, namely
[math]\displaystyle{
\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 \lt math\gt n }[/math] and [math]\displaystyle{ s }[/math] denote the normal and tangential directions on
[math]\displaystyle{ \partial \Delta }[/math] (where they exist) respectively and [math]\displaystyle{ \nu }[/math] is
Poisson's ratio.
Non-dimensional variables (denoted with an overbar) are introduced,
[math]\displaystyle{
(\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 \lt math\gt a }[/math] is a length parameter associated with the floe.
In non-dimensional variables, the equation for the ice floe
(plate_non) reduces to
[math]\displaystyle{ (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
\lt center\gt \lt math\gt
\beta = \frac{D}{g \rho a^4} \quad =and= \quad \gamma =
\frac{\rho_\Delta h}{ \rho a}.
The constants \lt math\gt \beta }[/math] and [math]\displaystyle{ \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)
[math]\displaystyle{ \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} &\lt \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]
and the Sommerfeld radiation condition
[math]\displaystyle{
\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 \lt math\gt \tilde{r}^2=x^2+y^2 }[/math] and [math]\displaystyle{ 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]\displaystyle{ d }[/math], the volume
occupied by the water changes, the vertical dimension being reduced to
[math]\displaystyle{ (-d,0) }[/math], (still denoted by [math]\displaystyle{ D }[/math]),
and the depth condition (water_depth) is replaced by the bed
condition,
[math]\displaystyle{
\frac{\partial \phi}{\partial z} = 0, \quad \mathbf{y} \in D,\: z=-d.
In water of finite depth, the positive real wavenumber \lt math\gt 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
[math]\displaystyle{
\kappa^* \tan \kappa^* d = - \frac{\alpha}{\beta \kappa^{*4} - \gamma
\alpha +1},
where \lt math\gt \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
[math]\displaystyle{
\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 \lt math\gt \kappa }[/math] smaller than [math]\displaystyle{ 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]\displaystyle{ 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,
[math]\displaystyle{ (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)
\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,
\lt math\gt \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,
[math]\displaystyle{ (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)
\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 \lt math\gt \nabla^4 }[/math], subject to the free edge boundary
conditions, is self-adjoint a thin plate must possess a set of modes [math]\displaystyle{ w^k }[/math]
which satisfy the free boundary conditions and the eigenvalue
equation
[math]\displaystyle{
\nabla^4 w^k = \lambda_k w^k.
The modes which correspond to different eigenvalues \lt math\gt \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.
[math]\displaystyle{
\int\limits_\Delta w^j (\mathbf{\xi}) w^k (\mathbf{\xi})
\d\sigma_{\mathbf{\xi}} = \delta _{jk},
where \lt math\gt \delta _{jk} }[/math] is the Kronecker delta. The eigenvalues [math]\displaystyle{ \lambda_k }[/math]
have the property that [math]\displaystyle{ \lambda_k \rightarrow \infty }[/math] as </math>k \rightarrow
\infty[math]\displaystyle{ 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 \lt math\gt \Delta }[/math].
We expand the displacement of the floe in a finite number of modes [math]\displaystyle{ M }[/math], i.e.
[math]\displaystyle{ (expansion)
w(\mathbf{x}) =\sum_{k=1}^{M} c_k w^k (\mathbf{x}).
\gt From the linearity of (int_eq_hs) the potential can be
written in the form
\lt center\gt \lt math\gt (expansionphi)
\phi(\mathbf{x}) =\phi^0(\mathbf{x}) + \sum_{k=1}^{M} c_k \phi^k (\mathbf{x}),
where \lt math\gt \phi^0 }[/math] and [math]\displaystyle{ \phi^k }[/math] respectively satisfy the integral equations
\begin{subequations} (phi)
[math]\displaystyle{ (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
\lt center\gt \lt math\gt (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) \d\sigma_{\mathbf{\xi}}.
\end{subequations}
The potential \lt math\gt \phi^0 }[/math] represents the potential due to the incoming wave
assuming that the displacement of the ice floe is zero. The potential
[math]\displaystyle{ \phi^k }[/math] represents the potential which is generated by the plate
vibrating with the [math]\displaystyle{ k }[/math]th mode in the absence of any input wave forcing.
We substitute equations (expansion) and (expansionphi) into
equation (plate_final) to obtain
[math]\displaystyle{ (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 \lt math\gt w^j }[/math] and integrate over
the plate (i.e. we take the inner product with respect to [math]\displaystyle{ w^j }[/math]) taking
into account the orthogonality of the modes [math]\displaystyle{ w^j }[/math] and obtain
[math]\displaystyle{ (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})
\d\sigma_{\mathbf{\xi}},
which is a matrix equation in \lt math\gt c_k }[/math].
Equation (final) cannot be solved without determining the modes of
vibration of the thin plate [math]\displaystyle{ w^k }[/math] (along with the associated
eigenvalues [math]\displaystyle{ \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]\displaystyle{ j }[/math]th floe is expanded in a
finite number of modes [math]\displaystyle{ M_j }[/math] (since the number of modes may not
necessarily be the same), i.e.
[math]\displaystyle{ (expansion_f)
w_j \left( \mathbf{x}\right) =\sum_{k=1}^{M_j} c_{jk} w_j^k (\mathbf{x}).
\gt From the linearity of (int_eq_hs) the potential can be
written in the form
\lt center\gt \lt math\gt (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 \lt math\gt \phi _{0} }[/math] and [math]\displaystyle{ \phi_j^k }[/math] respectively satisfy the integral equations
\begin{subequations} (phi_f)
[math]\displaystyle{ (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}) \d\sigma_{\mathbf{\xi}}
and
\lt center\gt \lt math\gt (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) \d\sigma_{\mathbf{\xi}}.
\end{subequations}
The potential \lt math\gt \phi_j^{0} }[/math] represents the potential due the incoming wave
assuming that the displacement of the ice floe is zero,
[math]\displaystyle{ \phi_j^k }[/math] represents the potential which is generated by the [math]\displaystyle{ j }[/math]th plate
vibrating with the [math]\displaystyle{ k }[/math]th mode in the absence of any input wave forcing. It
should be noted that [math]\displaystyle{ \phi_j^k(\mathbf{x}) }[/math] is, in general, non-zero
for [math]\displaystyle{ \mathbf{x}\in \Delta_{n} }[/math] (since the vibration of the [math]\displaystyle{ j }[/math]th
plate will result in potential under the [math]\displaystyle{ n }[/math]th plate).
We substitute equations (expansion_f) and (expansionphi_f) into
equation (plate_final) to obtain
[math]\displaystyle{ (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 \lt math\gt w_j^{l} }[/math] and integrate
over the plate (as before)
taking into account the orthogonality of the modes [math]\displaystyle{ w_j^l }[/math] and obtain
[math]\displaystyle{ (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}) \d\sigma_{\mathbf{\xi}},
which holds for all \lt math\gt j= 1, \ldots, N }[/math] and therefore gives a matrix
equation
