Kagemoto and Yue Interaction Theory

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.

We extend the finite depth interaction theory of kagemoto86 to water of infinite depth and bodies of arbitrary geometry. The sum over the discrete roots of the dispersion equation in the finite depth theory becomes an integral in the infinite depth theory. This means that the infinite dimensional diffraction transfer matrix in the finite depth theory must be replaced by an integral operator. In the numerical solution of the equations, this integral operator is approximated by a sum and a linear system of equations is obtained. We also show how the calculations of the diffraction transfer matrix for bodies of arbitrary geometry developed by goo90 can be extended to infinite depth, and how the diffraction transfer matrix for rotated bodies can be easily calculated. This interaction theory is applied to the wave forcing of multiple ice floes and a method to solve the full diffraction problem in this case is presented. Convergence studies comparing the interaction method with the full diffraction calculations and the finite and infinite depth interaction methods are carried out.

Introduction

The scattering of water waves by floating or submerged bodies is of wide practical importance. Although the problem is non-linear, if the wave amplitude is sufficiently small, the problem can be linearised. The linear problem is still the basis of the engineering design of most off-shore structures and is the standard model of geophysical phenomena such as the wave forcing of ice floes. While analytic solutions have been found for simplified problems (especially for simple geometries or in two dimensions) the full three-dimensional linear diffraction problem can only be solved by numerical methods involving the discretisation of the body's surface. The resulting linear system of equations has a dimension equal to the number of unknowns used in the discretisation of the body.

If more than one body is present, all bodies scatter the incoming waves. Therefore, the scattered wave from one body is incident upon all the others and, given that they are not too far apart, this notably changes the total incident wave upon them. Therefore, the diffraction calculation must be conducted for all bodies simultaneously. Since each body must be discretised this can lead to a very large number of unknowns. However, the scattered wavefield can be represented in an eigenfunction basis with a comparatively small number of unknowns. If we can express the problem in this basis, using what is known as an {\em Interaction Theory}, the number of unknowns can be much reduced, especially if there is a large number of bodies.

The first interaction theory that was not based on an approximation was the interaction theory developed by kagemoto86. Kagemoto and Yue found an exact algebraic method to solve the linear wave problem for vertically non-overlapping bodies in water of finite depth. The only restriction of their theory was that the smallest escribed circle for each body must not overlap any other body. The interaction of the bodies was accounted for by taking the scattered wave of each body to be the incident wave upon all other bodies (in addition to the ambient incident wave). Furthermore, since the cylindrical eigenfunction expansions are local, these were mapped from one body to another using Graf's addition theorem for Bessel functions. Doing this for all bodies, \citeauthor{kagemoto86} were able to solve for the coefficients of the scattered wavefields of all bodies simultaneously. The only difficulty with this method was that the solutions of the single diffraction problems had to be available in the cylindrical eigenfunction expansion of an outgoing wave. \citeauthor{kagemoto86} therefore only solved for axisymmetric bodies because the single diffraction solution for axisymmetric bodies was available in the required representation.

The extension of the Kagemoto and Yue scattering theory to bodies of arbitrary geometry was performed by goo90 who found a general way to solve the single diffraction problem in the required cylindrical eigenfunction representation. They used a representation of the finite depth free surface Green's function in the eigenfunction expansion of cylindrical outgoing waves centred at an arbitrary point of the water surface (above the body's mean centre position in this case). This Green's function was presented by black75 and further investigated by fenton78 who corrected some statements about the Green's function which Black had made. This Green's function is based on the cylindrical eigenfunction expansion of the finite depth free surface Green's function given by john2. The results of \citeauthor{goo90} were recently used by chakrabarti00 to solve for arrays of cylinders which can be divided into modules.

The development of the Kagemoto and Yue interaction theory was motivated by problems in off-shore engineering. However, the theory can also be applied to the geophysical problem of wave scattering by ice floes. At the interface of the open and frozen ocean an interfacial region known as the Marginal Ice Zone (MIZ) forms. The MIZ largely controls the interaction of the open and frozen ocean, especially the interaction through wave processors. This is because the MIZ consists of vast fields of ice floes whose size is comparable to the dominant wavelength, which means that it strongly scatters incoming waves. A method of solving for the wave response of a single ice floe of arbitrary geometry in water of infinite depth was presented by JGR02. The ice floe was modelled as a floating, flexible thin plate and its motion was expanded in the free plate modes of vibration. Converting the problem for the water into an integral equation and substituting the free modes, a system of equations for the coefficients in the modal expansion was obtained. However, to understand wave propagation and scattering in the MIZ we need to

understand the way in which large numbers of interacting

ice floes scatter waves. For this reason, we require an interaction theory. While the Kagemoto and Yue interaction theory could be used, their theory requires that the water depth is finite. While the water depth in the Marginal Ice Zone varies, it is generally located far from shore above the deep ocean. This means that the finite depth must be chosen large in order to be able to apply their theory. Furthermore, when ocean waves propagate beneath an ice floe the wavelength is increased so that it becomes more difficult to make the water depth sufficiently deep that it may be approximated as infinite. For this reason, in this paper we develop the equivalent interaction theory to Kagemoto and Yue's in infinite depth. Also, because of the complicated geometry of an ice floe, this interaction theory is for bodies of arbitrary geometry.

In the first part of this paper Kagemoto and Yue's interaction theory is extended to water of infinite depth. We represent the incident and scattered potentials in the cylindrical eigenfunction expansions and we use an analogous infinite depth Green's function to the one used by \citeauthor{goo90} \cite[given by][]{malte03}. We show how the infinite depth diffraction transfer matrices can be obtained with the use of this Green's function and we illustrate how the rotation of a body about its mean centre position in the plane can be accounted for without recalculating the diffraction transfer matrix.

In the second part of the paper, using \citeauthor{JGR02}'s single floe result, the full diffraction calculation for the motion and scattering from many interacting ice floes is calculated and presented. For two square interacting ice floes the convergence of the method obtained from the developed interaction theory is compared to the result of the full diffraction calculation. The solutions of more than two interacting ice floes and of other shapes in different arrangements are presented as well. We also compare the convergence of the finite depth and infinite depth methods in deep water.

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}, }[/math]

with discrete coefficients [math]\displaystyle{ 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, }[/math]

and the values of [math]\displaystyle{ 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. }[/math]

The incident potential upon body [math]\displaystyle{ \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}, }[/math]

with discrete coefficients [math]\displaystyle{ 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, }[/math]

where [math]\displaystyle{ {\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)}, }[/math]

given by black75 and fenton78, needs to be used instead of the infinite depth Green's function (green_inf). The elements of [math]\displaystyle{ {\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}) \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]

\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}\\ =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, }[/math]

with the density of the water [math]\displaystyle{ \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, }[/math]

where [math]\displaystyle{ 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}}, }[/math]

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

with

[math]\displaystyle{ \beta = \frac{D}{g \rho a^4} \quad =and= \quad \gamma = \frac{\rho_\Delta h}{ \rho a}. }[/math]

The constants [math]\displaystyle{ \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, }[/math]

\end{subequations} where [math]\displaystyle{ \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. }[/math]

In water of finite depth, the positive real wavenumber [math]\displaystyle{ 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}, }[/math]

where [math]\displaystyle{ \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}. }[/math]

For realistic values of the parameters, the effect of the plate is to make [math]\displaystyle{ \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) \mathrm{d}\sigma_\mathbf{\zeta}, \quad \mathbf{y} \in D. }[/math]

In the case of a shallow draft, the fact that the Green's function is symmetric and therefore satisfies the free surface boundary condition with respect to the second variable as well can be used to drastically simplify (int_eq). Due to the linearity of the problem the ambient incident potential can just be added to the equation to obtain the total water velocity potential, [math]\displaystyle{ \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) \mathrm{d}\sigma_\mathbf{\xi}, \quad \mathbf{x} \in \Delta. }[/math]

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

The modes which correspond to different eigenvalues [math]\displaystyle{ \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}) \mathrm{d}\sigma_{\mathbf{\xi}} = \delta _{jk}, }[/math]

where [math]\displaystyle{ \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}). }[/math]

>From the linearity of (int_eq_hs) the potential can be written in the form

[math]\displaystyle{ (expansionphi) \phi(\mathbf{x}) =\phi^0(\mathbf{x}) + \sum_{k=1}^{M} c_k \phi^k (\mathbf{x}), }[/math]

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

and

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

\end{subequations} The potential [math]\displaystyle{ \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. }[/math]

To solve equation (expanded) we multiply by [math]\displaystyle{ 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}) \mathrm{d}\sigma_{\mathbf{\xi}}, }[/math]

which is a matrix equation in [math]\displaystyle{ 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}). }[/math]

>From the linearity of (int_eq_hs) the potential can be written in the form

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

where [math]\displaystyle{ \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}) \mathrm{d}\sigma_{\mathbf{\xi}} }[/math]

and

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

\end{subequations} The potential [math]\displaystyle{ \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. }[/math]

To solve equation (expanded_f) we multiply by [math]\displaystyle{ 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}) \mathrm{d}\sigma_{\mathbf{\xi}}, }[/math]

which holds for all [math]\displaystyle{ j= 1, \ldots, N }[/math] and therefore gives a matrix equation for [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ \beta = 0.02 }[/math], mass [math]\displaystyle{ \gamma = 0.02 }[/math] and Poisson's ratio is chosen as [math]\displaystyle{ \nu=0.3333 }[/math]. The wavelength of the ambient incident wave is [math]\displaystyle{ \lambda = 2 }[/math]. Each ice floe has side length 2. The ambient wavefield is of unit amplitude and propagates in the [math]\displaystyle{ 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]\displaystyle{ O_1 }[/math] & [math]\displaystyle{ O_2 }[/math]\

[math]\displaystyle{ 3pt] 1 & \lt math\gt (-1.4,0) }[/math] & [math]\displaystyle{ (1.4,0) }[/math]\\

2 & [math]\displaystyle{ (0,-1.4) }[/math] & [math]\displaystyle{ (0,1.4) }[/math]\\ 3 & [math]\displaystyle{ (-1.4,-0.6) }[/math] & [math]\displaystyle{ (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}\[math]\displaystyle{ 0.4cm] \includegraphics[width=8.2cm]{int_n11_x0_y28_v5_b02_chi0}\\lt center\gt \lt math\gt 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. \lt center\gt \lt math\gt E_2 = \left( \, \int\limits_{\Delta} \big| w_{i}(\mathbf{x})-w_{f}(\mathbf{x}) \big|^2 \mathrm{d}\mathbf{x} \, \right)^{1/2}, }[/math]

where [math]\displaystyle{ w_{i} }[/math] and [math]\displaystyle{ w_{f} }[/math] are the solutions of the interaction method and the corresponding full diffraction calculation respectively. It would also be possible to compare other errors, the maximum difference of the solutions for example, but the results are very similar.

It is worth noting that the finite depth interaction method only converges up to a certain depth if used with the eigenfunction expansion of the finite depth Green's function (green_d). This is because of the factor [math]\displaystyle{ \alpha^2-k^2 }[/math] in the term of propagating modes of the Green's function. The Green's function can be rewritten by making use of the dispersion relation (eq_k) \cite[as suggested by][p. 26, for example]{linton01} and the depth restriction of the finite depth interaction method for bodies of arbitrary geometry can be circumvented.

The truncation parameters for the interaction methods will now be considered for both finite and infinite depth. The number of propagating modes and angular decaying components are free parameters in both methods. In finite depth, the number of decaying roots of the dispersion relation needs to be chosen while in infinite depth the discretisation of a continuous variable must be selected. In the infinite depth case we are free to choose the number of points as well as the points themselves. In water of finite depth, the depth can also be considered a free parameter as long as it is chosen large enough to account for deep water.

Truncating the infinite sums in the eigenfunction expansion of the outgoing water velocity potential for infinite depth with truncation parameters [math]\displaystyle{ T_H }[/math] and [math]\displaystyle{ T_K }[/math] and discretising the integration by defining a set of nodes, [math]\displaystyle{ 0\leq\eta_1 \lt \ldots \lt \eta_m \lt \ldots \lt \eta_{_{T_R}}\lt math\gt , with weights }[/math]h_m</math>, the potential for infinite depth can be approximated by

[math]\displaystyle{ \phi (r,\theta,z) &= \mathrm{e}^{\alpha z} \sum_{\nu = - T_H}^{T_H} A_{0\nu} H_\nu^{(1)} (\alpha r) \mathrm{e}^{\mathrm{i}\nu \theta}\\ &\quad + \sum_{m=1}^{T_R} h_m \, \psi (z,\eta_m) \sum_{\nu = - T_K}^{T_K} A_{\nu} (\eta_m) K_\nu (\eta_m r) \mathrm{e}^{\mathrm{i}\nu \theta}. }[/math]

In the following, the integration weights are chosen to be [math]\displaystyle{ h_m = 1/2\,(\eta_{m+1}-\eta_{m-1})\lt math\gt , }[/math]m=2, \ldots, T_R-1[math]\displaystyle{ and }[/math]h_1 = \eta_2-\eta_1[math]\displaystyle{ as well as }[/math]h_{_{T_R}} = \eta_{_{T_R}}-\eta_{_{T_R-1}}</math>, which corresponds to the mid-point quadrature rule. Different quadrature rules such as Gaussian quadrature could be considered. Although in general this would lead to better results, the mid-point rule allows a clever choice of the discretisation points so that the convergence with Gaussian quadrature is no better. In finite depth, the analogous truncation leads to

[math]\displaystyle{ \phi (r,\theta,z) &= \frac{\cosh k (z+d)}{\cosh kd} \sum_{\nu = - T_H}^{T_H} A_{0\nu} H_\nu^{(1)} (k r) \mathrm{e}^{\mathrm{i}\nu \theta}\\ & \quad + \sum_{m = 1}^{T_R} \frac{\cos k_m(z+d)}{\cos k_m d} \sum_{\nu = - T_K}^{T_K} A_{m\nu} K_\nu (k_m r) \mathrm{e}^{\mathrm{i}\nu \theta}. }[/math]

In both cases, the dimension of the diffraction transfer matrix, [math]\displaystyle{ \mathbf{B} }[/math], is given by [math]\displaystyle{ 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]\displaystyle{ T_H }[/math] to be 11 and [math]\displaystyle{ 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]\displaystyle{ T_H }[/math] and [math]\displaystyle{ T_K }[/math] are fixed (with the previously mentioned values) and [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ T_R }[/math] & discretisation of [math]\displaystyle{ \eta }[/math]\[math]\displaystyle{ 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 \lt math\gt 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]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ \lambda }[/math] has been chosen to be [math]\displaystyle{ 2 }[/math], the stiffness [math]\displaystyle{ \beta }[/math] and the mass [math]\displaystyle{ \gamma }[/math] of the ice floes to be 0.02 and Poisson's ratio [math]\displaystyle{ \nu }[/math] is [math]\displaystyle{ 0.3333 }[/math]. The ambient wavefield of amplitude 1 propagates in the positive direction of the [math]\displaystyle{ 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]\displaystyle{ 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}\[math]\displaystyle{ 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]