Difference between revisions of "Free-Surface Green Function for a Floating Elastic Plate"

From WikiWaves
Jump to navigationJump to search
 
(68 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 +
{{incomplete pages}}
 +
 
This is a special version of the free-surface Green function which applied when the [[Floating Elastic Plate]]
 
This is a special version of the free-surface Green function which applied when the [[Floating Elastic Plate]]
boundary condition applies at the free-surface
+
boundary condition applies at the free-surface. It reduced to the [[Free-Surface Green Function]] in the limit as the plate terms tend to zero and to the Green function for an infinite plate in the limit as the water terms tend to zero.
  
 
= Two Dimensions =
 
= Two Dimensions =
  
The Green function <math>G(x,z)</math> representing  
+
== Green function with singularity on the surface ==
outgoing waves as <math>|x|\rightarrow \infty</math> satisfies
+
For the case of a floating thin plate we almost always make the assumption of shallow draft and the Green function satisfies the following equations
 
+
The Green function <math>G(x,x^{\prime},z)</math> representing outgoing waves as <math>|x|\rightarrow \infty</math> satisfies
<math>
+
<center><math>
\nabla^2 G = 0, -h<z<0,  
+
\nabla^2 G = 0, -H<z<0,  
</math>
+
</math></center>
 
+
<center><math>
<math>
+
\frac{\partial G}{\partial z} =0, z=-H,
\frac{\partial G}{\partial z} =0, z=-h,
+
</math></center>
</math>
+
<center><math>
 
 
<math>
 
 
{\left( \beta \frac{\partial^4}{\partial x^4} -  
 
{\left( \beta \frac{\partial^4}{\partial x^4} -  
 
\gamma\alpha + 1\right)\frac{\partial G}{\partial z} - \alpha G}
 
\gamma\alpha + 1\right)\frac{\partial G}{\partial z} - \alpha G}
  = \delta(x), z=0,
+
  = \delta(x-x^{\prime}), z=0,
</math>
+
</math></center>
 
+
(note we are only considering singularities at the free surface). It is important to realise that these equations are based on the [[Frequency Domain Problem]] and that the exact form of the Green function is dependent on whether we have <math>e^{i\omega t}</math> or <math>e^{-i\omega t}</math> dependence. In what follows we assume <math>e^{i\omega t}</math> dependence.
(note we are only considering singularities at the free surface).  
 
 
 
This problem can be solved to give
 
 
 
<math>
 
G(x,z) = -i\sum_{n=-2}^\infty\frac{\sin{(k_n h)}\cos{(k(n)(z-h))}}{2\alpha C_n}e^{-k_n|x|},
 
</math>
 
  
 +
This problem can be solved (a solution is given below and also in [[Evans and Porter 2006]]) to give
 +
<center><math>
 +
G(x,x^{\prime},z) = -\sum_{n=-2}^\infty\frac{\sin{(k_n H)}\cos{(k_n(z+H))}}{2\alpha C(k_n)}e^{-k_n|x-x^{\prime}|},
 +
</math></center>
 
where
 
where
 
+
<center><math>
<math>
+
C(k_n)=\frac{1}{2}\left(h - \frac{(5\beta k_n^4 + 1 - \alpha\gamma)\sin^2{(k_n H)}}{\alpha}\right),
C_n=\frac{1}{2}\left(h + \frac{(5\beta k(n)^4 + 1 - \alpha\gamma)\sin^2{(k(n)h)}}{\alpha}\right),
+
</math></center>
</math>
 
 
 
 
and <math>k_n</math> are the solutions of the [[Dispersion Relation for a Floating Elastic Plate]],
 
and <math>k_n</math> are the solutions of the [[Dispersion Relation for a Floating Elastic Plate]],
 +
<center><math>
 +
\beta k^5 \sin(kH) - k \left(1 - \alpha \gamma \right) \sin(kH) =
 +
-\alpha \cos(kH)  \,
 +
</math></center>
 
with <math>n=-1,-2</math> corresponding to the complex solutions with positive real part,  
 
with <math>n=-1,-2</math> corresponding to the complex solutions with positive real part,  
<math>n=0</math> corresponding to the imaginary solution with positive imaginary part and
+
<math>n=0</math> corresponding to the imaginary solution with positive imaginary part, since we have assumed <math>e^{i\omega t}</math> dependence (it would be of negative imaginary part if we had chosen <math>e^{-i\omega t}</math> dependence)
 
<math>n>0</math> corresponding to the real solutions with positive real part.
 
<math>n>0</math> corresponding to the real solutions with positive real part.
 +
We can also write the Green function as
 +
<center><math>
 +
G\left(  r,z\right)  = -\omega\sum_{n=-2}^{\infty}
 +
\frac{R\left(  q_n\right)  }{q_n\sinh q_nH}\exp\left(  \mathrm{i}q_nr\right)
 +
\cosh q_n\left(  z+H\right)  .
 +
</math></center>
 +
where <math>r=|x-x^{\prime}|</math>, <math>q_n=ik_n</math> and we have used a different non-dimensionalisation so that
 +
the [[Dispersion Relation for a Floating Elastic Plate]] is written as
 +
<center><math>
 +
-q^5 \sinh(kH) - k \left(1 - m\omega^2 \right) \sinh(kH) =
 +
-\omega^2 \cosh(kH)  \,
 +
</math></center>
 +
and
 +
<center><math>
 +
R\left(  q_n\right)  =\frac{\omega^{2}q_n}{\omega^{2}\left(  5q_n^{4}+u\right)
 +
+H\left[  \left(  q_n^{5}+uq\right)  ^{2}-\omega^{4}\right]  }.
 +
</math></center>
  
= ThreeDimensions =
+
== Green function with singularity in the water ==
  
= Derivation of the Green function =
+
We can define the Green function <math>G(x,x^{\prime},z)</math>, representing
 
+
outgoing waves as <math>|x|\rightarrow \infty</math>, which satisfies
We present here a derivation of the Green function.
 
 
 
It uses the [[Mittag-Leffler Expansion for the  Floating Elastic Plate Dispersion Equation]]
 
 
 
In this chapter we consider the simple harmonic oscillation of an infinite
 
floating elastic plate using the method of solution reported in Fox
 
[[Fox00]], Fox and Chung [[report]], and Fox, Haskell and Chung
 
[[Fox01]]. We consider wave propagation in an infinite elastic plate that
 
is floating on water, i.e., the whole surface of the water is occupied by an
 
elastic plate, in our case an ice sheet satisfies the thin plate equation
 
((2-3)) given in chapter 2.
 
 
 
There are two main objectives of this chapter. First, we derive analytical
 
formulae for the deflection of an infinite ice sheet under a localized force.
 
Second, we scale the system of equations and the solutions to dimensionless
 
values so they are largely independent of physical parameter values within the ranges of our
 
interest. In order to derive the analytical solutions, we use the Fourier
 
transform in one and two dimensional space. Radially symmetric waves in a thin
 
plate without
 
the hydrodynamic support, have studied by Sneddon [[Sneddon45]] using the two
 
dimensional radial Fourier transform. This integral transform is commonly
 
called the Hankel transform
 
 
<center><math>
 
<center><math>
\hat{w}\left( \gamma\right)  =\frac{1}{2\pi}\int_{0}^{\infty}w\left(
+
\nabla^2 G = \delta(x-x^{\prime})\delta(z-z^{\prime}), -H<z<0,
r\right) rJ_{0}\left(  \gamma r\right)  dr.
 
 
</math></center>
 
</math></center>
A complete analytical solution for deflection of an infinite plate due to a
 
localized static load was derived by Wyman [[Wyman50]] by choosing the
 
bounded solutions of
 
 
<center><math>
 
<center><math>
w_{rr}+\frac{1}{r}w_{r}+\mathrm{i}w=0 (3-81)
+
\frac{\partial G}{\partial z} =0, z=-H,
 
</math></center>
 
</math></center>
whose solutions satisfy the radial form of the (non-dimensional) thin plate
 
equation
 
 
<center><math>
 
<center><math>
\left\frac{\partial^{2}}{\partial r^{2}}+\frac{1}{r}\frac{\partial
+
{\left( \beta \frac{\partial^4}{\partial x^4} -
}{\partial r}\right] ^{2}w+w=0.
+
\gamma\alpha + 1\right)\frac{\partial G}{\partial z} - \alpha G}
 +
  = 0, z=0.
 
</math></center>
 
</math></center>
 +
Note that this Green function reduces the [[Free-Surface Green Function]] in the limit as <math>\beta</math> and <math>\gamma</math> tend to zero.
  
The hydrodynamic effects that appear as <math>\phi_{t}</math> in the thin plate equation
+
This problem can be solved to give [[Evans and Porter 2003]]
force us to solve Laplace's equation, and thus the application of the Fourier
+
<center><math>
transform is not as straightforward as in the two examples cited (Sneddon
+
G(x,x^{\prime},z) = -\sum_{n=-2}^\infty\frac{\cos{(k_n(z^{\prime}+H))}\cos{(k_n(z+H))}}{2k_n C(k_n)}e^{-k_n|x-x^{\prime}|},
[[Sneddon45]] and Wyman [[Wyman50]]). Response of an infinite thin
+
</math></center>
elastic plate under fluid loading is studied in a series of papers by Crighton
+
where <math>
[[Crighton72, Crighton77, Crighton79]]. In the following sections, the
+
C(k_n)</math>
deflection of a plate is expressed by the infinite number of modes that exist
+
and <math>k_n</math> are as before.
in the vibrating plate floating on the water, instead of a finite number of
 
modes that represent the free oscillation of an elastic plate without
 
restraints. The fact that waves in an incompressible fluid can be expressed by
 
an infinite series of natural modes is shown by John [[fritz, fritz2]],
 
each mode being a Hankel function of the first kind. It is known that the
 
oscillation of an infinite plate can also be expressed in terms of Hankel
 
functions. Kheisin, in 1967 ([[Kheisin67]] chapter IV), studied the same
 
problem solved in this chapter. Kheisin derived the same dispersion equation
 
for the physical variables, and then studied the properties of the inverse
 
transform for the simple shallow water and static load cases, which are
 
approximated versions of the solutions reported here.
 
  
The most important consequence of being able to derive analytical solutions
+
== Relationship between the two Green functions ==
may be confirmation of the suitability of the scaling scheme for the range of
 
wave frequencies in which our geophysical interest lies (Fox [[Fox00]]).
 
Following Fox [[Fox00]], we show that
 
the characteristic length that appears in Wyman [[Wyman50]] for the
 
static-loading problem is also a natural length scale for a floating ice sheet
 
which is oscillating. In section (sec:2), we show that solutions scaled by
 
the characteristic length and the corresponding characteristic time are
 
independent of physical parameters, which enables us to find scaling laws
 
relating various scales of ice sheets.
 
  
In section (sec:1), we propose several methods of finding the
+
If we denote the Green function with the singularity at the free-surface by <math>G_s</math> and the Green function with
characteristic length, and thus the effective Young's modulus of a fast ice
+
the singularity in the water by <math>G_w</math>then we have the following relationship between them
sheet using the non-dimensional solutions. In Fox, et al. [[Fox01]],
+
<center><math>
characteristic length and effective Young's modulus are estimated from actual
+
G_s = -\frac{1}{\alpha}\left. \frac{\partial G_w}{\partial z^{\prime}} \right|_{z^{\prime} =0}
experimental data sets collected by Fox at\ McMurdo Sound, Antarctica.
+
</math></center>
  
 +
= Three Dimensions =
  
 +
<center><math>
 +
G\left(  r,z\right)  =-\frac{\omega}{2}\sum_{n=-2}^{\infty}
 +
\frac{R\left(  q_n\right)  }{\sinh q_nH}H_{0}^{\left(  1\right)  }\left(
 +
q_n r\right)  \cosh q_n\left(  z+H\right)
 +
</math></center>
  
===Fundamental solution===
+
= Derivation of the Green function =
  
In this section, we derive a fundamental solution for a localized forcing with
+
We present here a derivation of the Green function for both the two and three dimensional
a single frequency applied on the ice sheet that is occupying the whole
+
problem. The devivation uses the [[Mittag-Leffler Expansion for the Floating Elastic Plate Dispersion Equation]].
surface of the ocean. The fundamental solution for finite-depth water
+
This derivation is based on the PhD thesis of [[Hyuck Chung]].
was first derived by Fox (private communication) and then extended to the
 
infinite-depth case by the present author. We use the Fourier transform to
 
obtain the fundamental
 
solution which is expressed using an infinite series expansion over all wave
 
modes existing in a floating elastic plate. Then, we show the complete analytical
 
formulae of the coefficients of the modes and numerical computations of those
 
coefficients. The fundamental solution obtained here is normalized, that is,
 
the solution is expressed with dimensionless quantities by scaling the
 
distance and time. With our normalizing method, the BVP is simplified and
 
expressed without any physical values, which leads to the representation of
 
the solution being insensitive to certain physical values such as thickness of
 
the ice. In the following subsections, we will show details of derivation of
 
the spatial Fourier transform of <math>w\left(  x,y\right)  </math>, and then perform
 
analytical integrations of the inverse Fourier transform.
 
  
====Mathematical model====
+
==Previous Work ==
 +
 
 +
In the following sections, the
 +
deflection of a plate is expressed by the infinite number of modes that exist
 +
in the vibrating plate floating on the water, instead of a finite number of
 +
modes that represent the free oscillation of an elastic plate without
 +
restraints. The fact that waves in an incompressible fluid can be expressed by
 +
an infinite series of natural modes is shown by [[John 1949]] and [[John 1950]]
 +
and discussed in [[Cylindrical Eigenfunction Expansion]]
 +
Each mode being a Hankel function of the first kind in the horizontal
 +
(in the three-dimensional case) and a hyperbolic function
 +
in the vertical. [[Kheisin 1967]] chapter IV, studied the same
 +
problem solved in this chapter and derived the same dispersion equation
 +
for the physical variables. He studied the properties of the inverse
 +
transform only for the simple shallow water and static load cases, which are
 +
approximated versions of the full solutions which is presented here.
 +
 
 +
==Mathematical model==
  
 
In terms of non-dimensional variables the [[Floating Elastic Plate]] equation becomes
 
In terms of non-dimensional variables the [[Floating Elastic Plate]] equation becomes
 
<center><math>
 
<center><math>
 
\left[  \nabla^{4}-m\omega^{2}+1\right]  w+\mathrm{i}\omega
 
\left[  \nabla^{4}-m\omega^{2}+1\right]  w+\mathrm{i}\omega
\phi=p_{=a= }.  
+
\phi=p_n.  
 
</math></center>
 
</math></center>
Laplace's equation for the velocity potential and the bottom condition
+
where <math>p_n</math> is the pressure and <math>n</math> is the dimension. We assume
((2-12)) remain the same, i.e.,
+
that <math>p_2 = \delta(x) </math> and that <math>p_3 = \delta(x) \delta(y)</math>.
 +
<math>w</math> is the displacement of the surface which is related to the velocity
 +
potential by the equation
 +
<math>\frac{\partial\phi}{\partial z} = \frac{\partial w}{\partial t}</math>
 +
We also have Laplace's equation for the velocity potential and the bottom condition
 
<center><math>
 
<center><math>
 
\nabla^{2}\phi=0
 
\nabla^{2}\phi=0
Line 161: Line 144:
 
<math>\phi</math> for a given <math>\omega</math>.
 
<math>\phi</math> for a given <math>\omega</math>.
  
====Spatial Fourier transform====
+
==Spatial Fourier transform==
  
We solve the system Eqns.~((3-1)) to Eqns.~((3-3)) using the Fourier
+
We solve the system of equations using the [http://en.wikipedia.org/wiki/Fourier_Transform Fourier Transform] in <math>\left(  x,y\right)  </math>-plane for three-dimensions and with
transform in <math>\left(  x,y\right)  </math>-plane for point loading and on the
+
a [http://en.wikipedia.org/wiki/Fourier_Transform Fourier Transform] in
<math>x</math>-axis for line load. We choose the Fourier transform with respect to <math>x</math>
+
<math>x</math> for two dimensions. We choose the Fourier transform with respect to <math>x</math>
 
and <math>y</math> defined as
 
and <math>y</math> defined as
 
<center><math>
 
<center><math>
 
\hat{\phi}\left(  \alpha,k,z\right)  =\int_{-\infty}^{\infty}\int_{-\infty
 
\hat{\phi}\left(  \alpha,k,z\right)  =\int_{-\infty}^{\infty}\int_{-\infty
 
}^{\infty}\phi\left(  x,y,z\right)  e^{\mathrm{i}\left(  \alpha
 
}^{\infty}\phi\left(  x,y,z\right)  e^{\mathrm{i}\left(  \alpha
x+ky\right)  }dxdy (3-22)
+
x+ky\right)  }\mathrm{d}x\mathrm{d}y
 
</math></center>
 
</math></center>
and the inverse Fourier transform defined as
+
and the inverse [http://en.wikipedia.org/wiki/Fourier_Transform Fourier Transform] defined as
 
<center><math>
 
<center><math>
 
\phi\left(  x,y,z\right)  =\frac{1}{4\pi^{2}}\int_{-\infty}^{\infty}
 
\phi\left(  x,y,z\right)  =\frac{1}{4\pi^{2}}\int_{-\infty}^{\infty}
 
\int_{-\infty}^{\infty}\hat{\phi}\left(  \alpha,k,z\right)
 
\int_{-\infty}^{\infty}\hat{\phi}\left(  \alpha,k,z\right)
e^{-\mathrm{i}\left(  \alpha x+ky\right)  }d\alpha dk.
+
e^{-\mathrm{i}\left(  \alpha x+ky\right)  }\mathrm{d}\alpha \mathrm{d}k.
 
</math></center>
 
</math></center>
 
For the one-dimensional case, the definitions are
 
For the one-dimensional case, the definitions are
<center><math>\begin{matrix}
+
<center><math>
\hat{\phi}\left(  \alpha,z\right)   & =\int_{-\infty}^{\infty}\phi\left(
+
\hat{\phi}\left(  \alpha,z\right)   =\int_{-\infty}^{\infty}\phi\left(
x,z\right)  e^{\mathrm{i}\alpha x}dx, (3-36)\\
+
x,z\right)  e^{\mathrm{i}\alpha x}\mathrm{d}x,  
\phi\left(  x,z\right)   & =\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat{\phi
+
</math></center>
}\left(  \alpha,z\right)  e^{-\mathrm{i}\alpha x}d\alpha. (3-37)
+
<center><math>
\end{matrix}</math></center>
+
\phi\left(  x,z\right)   =\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat{\phi
We denote the spatial Fourier transform by using a hat over <math>w</math> and <math>\phi</math>.
+
}\left(  \alpha,z\right)  e^{-\mathrm{i}\alpha x}\mathrm{d}\alpha.  
 +
</math></center>
 +
We denote the spatial [http://en.wikipedia.org/wiki/Fourier_Transform Fourier Transform]
 +
by using a hat over <math>w</math> and <math>\phi</math>.
  
The Fourier transform of both sides of the Laplace's equation ((3-2))
+
The [http://en.wikipedia.org/wiki/Fourier_Transform Fourier Transform] of both sides of the Laplace's equation  
 
becomes an ordinary differential equation with respect to <math>z</math>,
 
becomes an ordinary differential equation with respect to <math>z</math>,
 
<center><math>
 
<center><math>
Line 202: Line 188:
 
function only of the magnitude of the Fourier transform variables,
 
function only of the magnitude of the Fourier transform variables,
 
<math>\gamma=\left\|  \mathbf{\gamma}\right\|  =\left|  \left|  \left(
 
<math>\gamma=\left\|  \mathbf{\gamma}\right\|  =\left|  \left|  \left(
\alpha,k\right)  \right|  \right|  <math>, hence we may now denote </math>\hat{\phi
+
\alpha,k\right)  \right|  \right|  </math>, hence we may now denote <math>\hat{\phi
}\left(  \alpha,k,z\right)  <math> by </math>\hat{\phi}\left(  \gamma,z\right)  <math>. We can reach the same conclusion using the fact that </math>w\left(  x,y\right)  <math> and </math>\phi\left(  x,y,z\right)  <math> are functions of </math>r=\sqrt{x^{2}+y^{2}}<math> thus the Fourier transforms must be functions of </math>\gamma=\sqrt{\alpha^{2}+k^{2}}</math>.
+
}\left(  \alpha,k,z\right)  </math> by <math>\hat{\phi}\left(  \gamma,z\right)  </math>. We can reach the same conclusion using the fact that <math>w\left(  x,y\right)  </math> and <math>\phi\left(  x,y,z\right)  </math> are functions of <math>r=\sqrt{x^{2}+y^{2}}</math> thus the Fourier transforms must be functions of <math>\gamma=\sqrt{\alpha^{2}+k^{2}}</math>.
  
 
We find the unknown coefficients <math>A</math> and <math>B</math> from the Fourier transformed ocean floor condition that <math>\left.  \hat{\phi}_{z}\right|  _{z=-H}=0</math> to be <math>A\left(  \gamma\right)  =Ce^{\gamma H}</math>, and <math>B\left(  \gamma\right) =Ce^{-\gamma H}</math>. Thus, we obtain the depth dependence of the Fourier transform of the potential
 
We find the unknown coefficients <math>A</math> and <math>B</math> from the Fourier transformed ocean floor condition that <math>\left.  \hat{\phi}_{z}\right|  _{z=-H}=0</math> to be <math>A\left(  \gamma\right)  =Ce^{\gamma H}</math>, and <math>B\left(  \gamma\right) =Ce^{-\gamma H}</math>. Thus, we obtain the depth dependence of the Fourier transform of the potential
 
<center><math>
 
<center><math>
 
\hat{\phi}\left(  \gamma,z\right)  =\hat{\phi}\left(  \gamma,0\right)
 
\hat{\phi}\left(  \gamma,z\right)  =\hat{\phi}\left(  \gamma,0\right)
\frac{\cosh\gamma\left(  z+H\right)  }{\cosh\gamma H}. (3-5)
+
\frac{\cosh\gamma\left(  z+H\right)  }{\cosh\gamma H}.  
 
</math></center>
 
</math></center>
At the surface, <math>z=0</math>, differentiating both sides of Eqn.~((3-5)) with respect to <math>z</math>,\ the vertical component of the velocity is
+
At the surface, <math>z=0</math>, differentiating both sides with respect to <math>z</math> the vertical component of the velocity is
 
<center><math>
 
<center><math>
 
\hat{\phi}_{z}\left(  \gamma,0\right)  =\hat{\phi}\left(  \gamma,0\right)
 
\hat{\phi}_{z}\left(  \gamma,0\right)  =\hat{\phi}\left(  \gamma,0\right)
\gamma\tanh\gamma H. (3-27)
+
\gamma\tanh\gamma H.  
 
</math></center>
 
</math></center>
 
Using this relationship to substitute for <math>\hat{\phi}_{z}</math> in the non-dimensional form of the kinematic condition, we find that
 
Using this relationship to substitute for <math>\hat{\phi}_{z}</math> in the non-dimensional form of the kinematic condition, we find that
 
<center><math>
 
<center><math>
 
\hat{\phi}\left(  \gamma,0\right)  =\frac{\mathrm{i}\omega\hat
 
\hat{\phi}\left(  \gamma,0\right)  =\frac{\mathrm{i}\omega\hat
{w}\left(  \gamma\right)  }{\gamma\tanh\gamma H}. (3-6)
+
{w}\left(  \gamma\right)  }{\gamma\tanh\gamma H}.  
 
</math></center>
 
</math></center>
 
Thus, once we find <math>\hat{w}\left(  \gamma\right)  </math> then we can find
 
Thus, once we find <math>\hat{w}\left(  \gamma\right)  </math> then we can find
<math>\hat{\phi}\left(  \gamma,z\right)  </math> using Eqn.~((3-27)) and
+
<math>\hat{\phi}\left(  \gamma,z\right)  </math>.
Eqn.~((3-6)).
 
  
The Fourier transform of the plate equation ((3-1)) is also an algebraic
+
The Fourier transform of the plate equation is also an algebraic
 
equation in the parameter <math>\gamma</math>
 
equation in the parameter <math>\gamma</math>
 
<center><math>
 
<center><math>
Line 230: Line 215:
 
\hat{\phi}=1.
 
\hat{\phi}=1.
 
</math></center>
 
</math></center>
Hence, substituting Eqn.~((3-6)), we have the single algebraic equation
+
Hence, we have the single algebraic equation
 
for each spatial Fourier component of <math>w</math>
 
for each spatial Fourier component of <math>w</math>
 
<center><math>
 
<center><math>
 
\left(  \gamma^{4}+1-m\omega^{2}-\frac{\omega^{2}}{\gamma\tanh\gamma
 
\left(  \gamma^{4}+1-m\omega^{2}-\frac{\omega^{2}}{\gamma\tanh\gamma
H}\right)  \hat{w}=1. (3-7)
+
H}\right)  \hat{w}=1.  
 
</math></center>
 
</math></center>
 
Notice that we have used the fact that the Fourier transform of the delta
 
Notice that we have used the fact that the Fourier transform of the delta
Line 243: Line 228:
 
<center><math>
 
<center><math>
 
\hat{w}\left(  \gamma\right)  =\frac{1}{d\left(  \gamma,\omega\right)
 
\hat{w}\left(  \gamma\right)  =\frac{1}{d\left(  \gamma,\omega\right)
} (3-8)
+
}  
 
</math></center>
 
</math></center>
 
where
 
where
 
<center><math>
 
<center><math>
 
d\left(  \gamma,\omega\right)  =\gamma^{4}+1-m\omega^{2}-\frac{\omega^{2}
 
d\left(  \gamma,\omega\right)  =\gamma^{4}+1-m\omega^{2}-\frac{\omega^{2}
}{\gamma\tanh\gamma H}. (3-9)
+
}{\gamma\tanh\gamma H}.  
 
</math></center>
 
</math></center>
We call the function <math>d\left(  \gamma,\omega\right)  </math>\ the dispersion
+
Where <math>d\left(  \gamma,\omega\right)  </math> in the  
function, and the associated algebraic equation <math>d\left(  \gamma
+
[[Dispersion Relation for a Floating Elastic Plate]] function, and the associated algebraic equation <math>d\left(  \gamma
 
,\omega\right)  =0</math> the dispersion equation for waves propagating through an
 
,\omega\right)  =0</math> the dispersion equation for waves propagating through an
ice sheet. This dispersion equation (and the Fourier
+
ice sheet. This dispersion relation (and the Fourier
transform Eqn.~((3-8))) has been previously derived by Kheisin
+
transform of it) was derived by  
([[Kheisin67]] chapter IV), Fox and Squire [[fox3]].
+
[[Kheisin 1967]] chapter IV, and [[Fox and Squire 1994]].
 
 
\begin{figure}[tbh]
 
\begin{center}
 
\includegraphics[height=4.3405cm,width=13.6476cm]{rootfind.eps}
 
\caption{Illustrates how to find the positions of the roots of the dispersion equation. On the right, the real roots <math>\pm q_{=T= }</math> and on the left the imaginary roots <math>iq_{1},iq_{2},iq_{3},....</math>. Dotted curves are <math>\tanh</math> function on the left and <math>\tan</math> function on the right. Solid curves are the curves of Eqn. ((3-38)).}
 
(fig:3-12)
 
\end{center}
 
\end{figure}
 
  
Our task now is to derive the inverse Fourier transform of Eqn.~((3-8)). We notice that the roots of the dispersion equation for a fixed <math>\omega</math>, are the poles of the function <math>\hat{w}\left(  \gamma\right)  </math>, which are necessary for calculation of integrals involved in the inverse Fourier transform. The roots of the dispersion equation are computed numerically using computer codes in MatLab given by Fox and Chung [[report]]. The root-finder computer code is due to Fox and was initially written to implement the mode matching technique in Fox and Squire [[fox3]]. Later, we will use the observation that the dispersion equation <math>d\left(  \gamma,\omega\right)  </math> is an even function of <math>\gamma</math>, and hence if <math>q</math> is a root then so is <math>-q</math>. Fig.~((fig:3-12)) shows how the real and imaginary roots can be found from the curves of the <math>\tan</math> and <math>\tanh</math> functions intersecting the curves
+
==Inverse Fourier Transform ==
<center><math>
 
\frac{\pm\omega^{2}}{\gamma\left(  \gamma^{4}+1-m\omega^{2}\right)
 
}. (3-38)
 
</math></center>
 
  
We notice that from observing the curves in Fig.~((fig:3-12)), there are a pair of two real roots <math>\pm q_{=T= }</math> (<math>q_{=T= }>0</math>) and an infinite number of pure imaginary roots <math>\left\{ \pm\mathrm{i}q_{n}\right\} _{n=1,2,...}</math>, (<math>q_{n}>0</math>). We notice that for typical cases <math>m\omega^{2} \leq1</math>, an imaginary root <math>\mathrm{i}q_{n}</math> is found in each interval <math>\left(  n-\frac{1}{2}\right)  \pi<\gamma_{n}H<n\pi</math> as shown in Fig.~((fig:3-12)). For fixed <math>\omega\neq0</math>, it is known (Fox and Squire [[fox3]], Chung and Fox [[chung]]) that four complex roots occur as plus and minus complex-conjugate pairs <math>\pm q_{=D= }</math> and <math>\pm q_{=D=  }^{}</math> (<math>Re \left(  q_{=D= }\right)  >0</math> and <math>Im \left(  q_{=D= }\right)  >0</math>).
+
Our task now is to derive the inverse Fourier transform). We notice that the roots of the
 +
[[Dispersion Relation for a Floating Elastic Plate]] for a fixed <math>\omega</math>, are the poles of the function <math>\hat{w}\left(  \gamma\right)  </math>, which are necessary for calculation of integrals involved in the inverse Fourier transform. There are a pair of two real roots <math>\pm q_0</math>, an infinite number of pure imaginary roots  
 +
<math>  \pm q_{n},\,n\in\mathbb{N}</math>, (<math>\mathrm{Im} q_{n}>0</math>) plus four complex roots  
 +
(in all physical situations) occur at plus and minus complex-conjugate pairs <math>\pm q_{-1}</math> and  
 +
<math>\pm q_{-2}</math> with <math>\mathrm{Im} \left(  q_{-1,-2}\right)  >0</math>).
  
 
We note that the dispersion equation represents a relationship between the spatial wavenumbers <math>\gamma</math> and the radial frequency <math>\omega</math>, which is how the name "dispersion" came about.
 
We note that the dispersion equation represents a relationship between the spatial wavenumbers <math>\gamma</math> and the radial frequency <math>\omega</math>, which is how the name "dispersion" came about.
Line 282: Line 259:
 
which is also a function of <math>\gamma</math> only and has the same poles as <math>\hat{w}\left(  \gamma\right)  </math> does.
 
which is also a function of <math>\gamma</math> only and has the same poles as <math>\hat{w}\left(  \gamma\right)  </math> does.
  
The Fourier transform of the differential equations can be interpreted as wave-like forcing of the ice sheet, i.e.,
+
We calculate the displacement <math>w\left(  x,y\right)  </math> using the
<center><math>
+
inverse Fourier transform of <math>\hat{w}\left(  \gamma\right)  </math>. Since <math>\hat{w}\left(  \gamma\right)  </math> is radially symmetric, the inverse transform may be written ([[Bracewell 1965]])
p_{=a= }\left(  x,y,t\right)  =\hat{p}_{=a= }\left(  \alpha
 
,k,\omega\right)  e^{\mathrm{i}\left(  \omega t-\alpha x-ky\right)  }
 
</math></center>
 
where <math>\hat{p}_{=a= }\left(  \alpha,k,\omega\right)  </math> is the amplitude of the wave-like forcing or a Fourier transform of the wave-like force which we set to be one. Therefore, superposition or the inverse Fourier transform of the solution under this wave-like force represents the response of a plate to a localized force expressed by the delta function in Eqn.~((3-23)). Since the waves are radially symmetric, <math>p_{=a= }</math> can be written as a function of <math>r</math>,
 
<center><math>
 
p_{=a= }\left(  r,t\right)  =\hat{p}_{=a= }\left(  \gamma
 
,\omega\right)  J_{0}\left(  \gamma r\right)  e^{\mathrm{i}\omega t}
 
</math></center>
 
where <math>\hat{p}_{=a= }\left(  \gamma,\omega\right)  </math> is the Fourier component of radially symmetric wave-like function.
 
 
 
====The inverse Fourier transform (sec:3)====
 
 
 
We calculate the displacement <math>w\left(  x,y\right)  </math> by performing the two dimensional inverse Fourier transform of <math>\hat{w}\left(  \gamma\right)  </math> in Eqn.~((3-8)) and, since <math>\hat{w}\left(  \gamma\right)  </math> is radially symmetric, the inverse transform may be written (Bracewell [[Bracewell65]])
 
 
<center><math>
 
<center><math>
w_{=P= }\left(  r\right)  =\frac{1}{2\pi}\int_{0}^{\infty}\hat{w}\left(
+
w\left(  r\right)  =\frac{1}{2\pi}\int_{0}^{\infty}\hat{w}\left(
\gamma\right)  \gamma J_{0}\left(  \gamma r\right)  d\gamma (3-11)
+
\gamma\right)  \gamma J_{0}\left(  \gamma r\right)  \mathrm{d}\gamma  
 
</math></center>
 
</math></center>
 +
for the three dimensional problem,
 
where <math>J_{0}</math> is Bessel function and <math>r</math> is the distance from the point of
 
where <math>J_{0}</math> is Bessel function and <math>r</math> is the distance from the point of
forcing. The response to line forcing is given by the inverse Fourier
+
forcing. For the two-dimensional problem the inverse Fourier
 
transform of <math>\hat{w}\left(  \gamma\right)  </math> in the <math>x</math>-axis and since
 
transform of <math>\hat{w}\left(  \gamma\right)  </math> in the <math>x</math>-axis and since
 
<math>\hat{w}\left(  \gamma\right)  </math> is an even function, this is
 
<math>\hat{w}\left(  \gamma\right)  </math> is an even function, this is
 
<center><math>
 
<center><math>
w_{=L= }\left(  r\right)  =\frac{1}{\pi}\int_{0}^{\infty}\hat{w}\left(
+
w\left(  r\right)  =\frac{1}{\pi}\int_{0}^{\infty}\hat{w}\left(
\gamma\right)  \cos\left(  \gamma r\right)  d\gamma (3-12)
+
\gamma\right)  \cos\left(  \gamma r\right)  \mathrm{d}\gamma  
 
</math></center>
 
</math></center>
where again <math>r=\left|  x\right|  </math> is the distance from the line of forcing. Note that the factors <math>1/2\pi</math> and <math>1/\pi</math> result from the form of the Fourier transform in two and one dimensional spaces that we defined in Eqn.~((3-22)) and Eqn.~((3-36)).
+
where again <math>r=\left|  x\right|  </math>. Note that the factors <math>1/2\pi</math> and <math>1/\pi</math> result from the form of the Fourier transform in two and one dimensional spaces.
  
The integrals of Eqn.~((3-11)) and Eqn.~((3-12)) are calculated using the singularities of <math>\hat{w}</math>, i.e., the roots of the dispersion equation. First, since <math>\hat{w}\left(  \gamma\right)  </math> is an even fractional function that equals zero when <math>\gamma=0</math> and is bounded in the whole plane except in regions around its poles, we find that <math>\hat{w}\left(  \gamma\right)  </math> can be expressed using the Mittag-Leffler expansion in section 3.2. Pairing each pole <math>q</math> with its negative counterpart <math>-q</math>, gives
+
The integrals are calculated using the [[Mittag-Leffler Expansion for the Floating Elastic Plate Dispersion Relation]]
 +
where we show that
 
<center><math>
 
<center><math>
\hat{w}\left(  \gamma\right)  =\sum_{q\in K^{}}\frac{2qR\left(
+
\hat{w}\left(  \gamma\right)  =\sum_{n=-2}^{\infty}\frac{2q_nR\left(
q\right)  }{\gamma^{2}-q^{2}} (3-50)
+
q_n\right)  }{\gamma^{2}-q_n^{2}}  
 
</math></center>
 
</math></center>
where <math>R\left(  q\right)  </math> is the residue of <math>\hat{w}\left(  \gamma\right)  </math> at <math>\gamma=q</math>. We denoted the set of roots of the dispersion equation with positive imaginary part and a positive real root by <math>K^{}</math>. Thus, this set is <math>K^{}=\left\{  q_{=T= },q_{=D= },-q_{=D=  }^{},\mathrm{i}q_{1},\mathrm{i}q_{2},\mathrm{i} q_{3},\cdots\right\}  </math>. Note that the rest of the roots of the dispersion equation are the negative of the values in <math>K^{}</math>. By substituting this expansion into the integrals in Eqn.~((3-11)) and Eqn.~((3-12)), we are able to perform the integration and write each result as a summation over the roots <math>q\in K^{}</math>.
+
where <math>R\left(  q_n\right)  </math> is the residue of <math>\hat{w}\left(  \gamma\right)  </math> at <math>\gamma=q_n</math>
 
+
given by
\begin{figure}[tbh]
 
\begin{center}
 
\includegraphics[height=5.0918cm,width=7.5037cm]{polyres.eps}
 
\caption{Graphs of <math>\left|  R\left(  \mathrm{i}q\right)  \right|  </math> as function of real number <math>q</math>, with <math>\tanh</math>-functions (solid line) and polynomial expression (dashed line).}
 
(fig:3-31)
 
\end{center}
 
\end{figure}
 
 
 
The residues <math>R\left(  q\right)  </math> can be calculated using the usual formula.
 
Since each of the poles of <math>\hat{w}\left(  \gamma\right)  </math> is simple, the
 
residue <math>R\left(  q\right)  </math> at a pole <math>q</math> can be found using the expression
 
<center><math>\begin{matrix}
 
R\left(  q\right)  & =\left[  \left.  \frac{d}{d\gamma}d\left(  \gamma
 
,\omega\right)  \right|  _{\gamma=q}\right]  ^{-1}\\
 
& =\left[  4q^{3}+\omega^{2}\left(  \frac{qH+\tanh qH-qH\tanh^{2}qH}
 
{q^{2}\tanh^{2}qH}\right)  \right]  ^{-1}. (3-16)
 
\end{matrix}</math></center>
 
As each pole <math>q</math> is a root of the dispersion equation, we may substitute
 
<math>\tanh qH=\omega^{2}/\left(  q^{5}+uq\right)  </math>, where for brevity we have
 
defined <math>u=\left(  1-m\omega^{2}\right)  </math>. The residue may then be given as
 
the rational function of the pole
 
 
<center><math>
 
<center><math>
R\left(  q\right)  =\frac{\omega^{2}q}{\omega^{2}\left(  5q^{4}+u\right)
+
R\left(  q_n\right)  =\frac{\omega^{2}q_n}{\omega^{2}\left(  5q_n^{4}+u\right)
+H\left[  \left(  q^{5}+uq\right)  ^{2}-\omega^{4}\right]  }. (3-17)
+
+H\left[  \left(  q_n^{5}+uq_n\right)  ^{2}-\omega^{4}\right]  }.
</math></center>
+
</math></center>.
This form avoids calculation of the hyperbolic tangent which becomes small at the imaginary roots. This causes numerical round-off problems and rapid variation in computing Eqn.~((3-16)) since <math>q_{n}H</math> tends to <math>n\pi</math> as <math>n</math> becomes large, which makes <math>\tan q_{n}H</math> become small. Fig.~((fig:3-31)) shows the graphs of the two expressions, Eqn.~((3-16)) and Eqn.~((3-17)) of the residue for imaginary argument. The imaginary roots <math>\mathrm{i}q_{1},\mathrm{i}q_{2},\mathrm{i}q_{3},...</math> are located where the two curves in Fig.~((fig:3-31)) coincide (the spiky parts of the solid curve). There are two such points at a spike and the root is the one on the left. Eqn.~((3-16)), from the direct calculation, is a rapidly varying function near the roots <math>\left\{  \mathrm{i} q_{n}\right\}  _{n=1}^{\infty}</math>, hence a small numerical error in the values of the roots will result in a large error in the residue. In contrast, Eqn.~((3-17)) gives us a smooth function, and the resulting calculation of the residue is stable.
 
  
Using the identities (Abramowitz and Stegun [[Stegun]] formula 11.4.44 with <math>\nu=0</math>, <math>\mu=0</math>, <math>z=-\mathrm{i}q</math> and <math>a=r</math>)
+
Using the identities [[Abramowitz and Stegun 1964]] formula 11.4.44 with <math>\nu=0</math>, <math>\mu=0</math>, <math>z=-\mathrm{i}q</math> and <math>a=r</math>)
 
<center><math>
 
<center><math>
 
\int_{0}^{\infty}\frac{\gamma}{\gamma^{2}-q^{2}}J_{0}\left(  \gamma r\right)
 
\int_{0}^{\infty}\frac{\gamma}{\gamma^{2}-q^{2}}J_{0}\left(  \gamma r\right)
d\gamma=K_{0}\left(  -\mathrm{i}qr\right)  (3-24)
+
\mathrm{d}\gamma=K_{0}\left(  -\mathrm{i}qr\right)   
</math></center>
 
for <math>Im q>0</math>, <math>r>0</math>, where <math>K_{0}</math> is a modified Bessel function and an identity between the modified Bessel function, we notice that <math>q_{=T= }</math> term of Eqn.~((3-50)) may pose a problem. However, considering that <math>q_{=T= }=\lim_{\varepsilon\searrow0}\left(  q_{=T=  }+\mathrm{i}\varepsilon\right)  </math>, we are able to apply Eqn.~((3-24)) to all terms of Eqn.~((3-50)). Alternatively, we may put an additional imaginary term <math>\mathrm{i}\beta\omega</math>, <math>\beta>0</math>, representing damping, to the dispersion equation
 
<center><math>
 
\gamma^{4}+\mathrm{i}\beta\omega+1-m\omega^{2}-\frac{\omega^{2}
 
}{\gamma\tanh\gamma H}.
 
 
</math></center>
 
</math></center>
Every element in <math>K^{}</math> has positive imaginary part. Addition of the damping term proportional to <math>w_{t}</math> ensures that the inverse transform satisfies the radiation condition. We omit the damping term here for the sake of algebraic simplicity.
+
for <math>Im q>0</math>, <math>r>0</math>, where <math>K_{0}</math> is a modified Bessel function and an identity between the modified Bessel function.
  
We can calculate the integral of the inverse Fourier transform of <math>\hat{w}</math> and, using the identity between <math>K_{0}</math> and Hankel function of the first kind (Abramowitz and Stegun [[Stegun]] formula 9.6.4), <math>K_{0}\left( \zeta\right)  =\frac{\mathrm{i}\pi}{2}H_{0}^{\left(  1\right)  }\left( \mathrm{i}\zeta\right)  </math> for <math>Re \zeta\geq0</math>, the displacement for point forcing may be written in the equivalent forms
+
We can calculate the integral of the inverse Fourier transform of <math>\hat{w}</math> and, using the identity between <math>K_{0}</math> and Hankel function of the first kind ([[Abramowitz and Stegun 1964]] formula 9.6.4),  
<center><math>\begin{matrix}
 
w_{=P= }\left(  r\right)  & =\frac{1}{\pi}\sum_{q\in K^{}
 
}qR\left(  q\right)  K_{0}\left(  -\mathrm{i}qr\right)
 
(3-13)\\
 
& =\frac{\mathrm{i}}{2}\sum_{q\in K^{}}qR\left(  q\right)
 
H_{0}^{\left(  1\right)  }\left(  qr\right)  (3-14)
 
\end{matrix}</math></center>
 
where the subscript P of <math>w_{=P= }</math> indicates the response to a point load.
 
 
 
The identity (Erdelyi [[bateman]] formula 1.2 (11) with <math>x=\gamma</math>, <math>y=r</math>, and <math>a=-\mathrm{i}q</math>)
 
 
<center><math>
 
<center><math>
\int_{0}^{\infty}\frac{\cos\left( \gamma r\right)  }{\gamma^{2}-q^{2}}
+
K_{0}\left( \zeta\right)  =\frac{\mathrm{i}\pi}{2}H_{0}^{\left(  1\right)  }\left( \mathrm{i}\zeta\right)   
d\gamma=-\frac{\pi}{\mathrm{i}2q}\exp\left( \mathrm{i}qr\right)
 
  (3-25)
 
 
</math></center>
 
</math></center>
for <math>Im q>0</math>, <math>r>0</math>. The above integration gives the surface displacement for line forcing
+
for <math>Re \zeta\geq0</math>, the displacement for three-dimensions may be written in the equivalent forms
 
<center><math>
 
<center><math>
w_{=L= }\left(  r\right) =\mathrm{i}\sum_{q\in K^{}
+
w\left(  r\right)   =\frac{1}{\pi}\sum_{n=-2}^{\infty}
}R\left(  q\right)  \exp\left(  \mathrm{i}qr\right)  (3-15)
+
q_n R\left(  q_n\right)  K_{0}\left(  -\mathrm{i}q_n r\right)
 
</math></center>
 
</math></center>
where the subscript L of <math>w_{=L= }</math> indicates the response to a line load.
 
 
====Modal expansion of the solutions====
 
 
We have shown that the fundamental solution for finite water depth <math>H</math> is
 
found by first finding the roots in the upper-half plane, <math>K^{}</math>,
 
of the dispersion Eqn.~((3-9)). Numerical calculation is achieved by
 
truncating the sum after some finite number of
 
roots. Straightforward computer code (in MatLab) to find the roots has been
 
given by Fox and Chung [[report]]. After calculating the residue for each
 
root given in Eqn.~((3-17)), the sum in Eqn.~((3-14)) can be
 
calculated by separating it into
 
 
<center><math>
 
<center><math>
w_{=P= }\left(  r\right)  =\frac{\mathrm{i}}{2}q_{=T=
+
=\frac{\mathrm{i}}{2}\sum_{n=-2}^{\infty}q_n R\left(  q_n\right)
}R_{=T= }H_{0}^{\left(  1\right)  }\left(  q_{=T= }r\right)
+
H_{0}^{\left(  1\right)  }\left(  q_n r\right).
-Im \left[  q_{=D= }R_{=D= }H_{0}^{\left(  1\right)
 
}\left(  q_{=D= }r\right)  \right] +\frac{1}{\pi}\sum_{n=1}^{\infty
 
}\mathrm{i}q_{n}R_{n}K_{0}\left(  q_{n}r\right) (3-18)
 
 
</math></center>
 
</math></center>
We used the identities <math>-\mathrm{i}\left(  -q_{=D= }\right)
 
=\left(  -\mathrm{i}q_{=D= }\right)  ^{}<math>, </math>R\left(  q^{
 
}\right)  =\left(  R\left(  q\right)  \right)  ^{}<math> and </math>H_{0}^{\left(
 
1\right)  }\left(  -\zeta^{}\right)  =-\left(  H_{0}^{\left(  1\right)
 
}\left(  \zeta\right)  \right)  ^{}</math>, and denote the residues for the
 
poles in <math>K^{}</math> by <math>R_{=T= }=R\left(  q_{=T= }\right)  </math>,
 
<math>R_{=D= }=R\left(  q_{=D= }\right)  </math>, and <math>R_{n}=R\left(
 
\mathrm{i}q_{n}\right)  <math>, </math>n\in\mathbf{N}</math>, respectively. Note that
 
<math>R\left(  -\gamma_{=D= }^{}\right)  \rightarrow-R_{=D= }^{}</math>
 
as <math>\beta\rightarrow0</math>. Since, from Eqn.~((3-17)), <math>\left|  qR\left(
 
q\right)  \right|  \propto\left|  q\right|  ^{-8}<math> as </math>\left|  q\right|  </math>
 
increases, the sum may be terminated after a relatively small number of terms
 
without significant error, as shown in Fig.~((fig:rsd)). We also notice
 
that the smaller the depth <math>H</math>, the fewer number of roots we need to achieve
 
the desired accuracy, since for a given <math>n</math>, <math>q_{n}</math> increases as <math>H</math> becomes small.
 
  
\begin{figure}[tbh]
+
In two-dimensions the identity ([[Erdelyi et. al.1954]] formula 1.2 (11) with <math>x=\gamma</math>, <math>y=r</math>, and <math>a=-\mathrm{i}q</math>)
\begin{center}
 
\includegraphics[height=7.15cm,width=7.3126cm]{rsd.eps}
 
\caption{Log-log plot of the residues corresponding to the evanescent modes. +, * and <math>\diamondsuit</math> indicate for <math>\omega=0.1,1.0,10</math> respectively. The water depth is <math>H=20\pi</math> (deep water).}
 
(fig:rsd)
 
\end{center}
 
\end{figure}
 
 
 
The surface displacement for line forcing is given by the sum
 
 
<center><math>
 
<center><math>
w_{=L= }\left(  r\right)  =\mathrm{i}R_{=T= }\exp\left(
+
\int_{0}^{\infty}\frac{\cos\left(  \gamma r\right)  }{\gamma^{2}-q^{2}}
\mathrm{i}q_{=T= }r\right)  -2Im \left[  R_{=D=
+
\mathrm{d}\gamma=-\frac{\pi}{\mathrm{i}2q}\exp\left(  \mathrm{i}qr\right)
}\exp\left(  \mathrm{i}q_{=D= }r\right)  \right]
 
+\mathrm{i}\sum_{n=1}^{\infty}R_{n}\exp\left(  -q_{n}r\right)
 
. (3-19)
 
</math></center>
 
We have written the solutions ((3-18)) and ((3-19)) in terms of the
 
travelling, damped travelling and evanescent mode in order to emphasize the
 
behavior of each mode.
 
 
 
We consider here the energy flux due to the wave produced by the force at
 
<math>r=0</math>. From the equation of motion ((2-41)) and <math>\mathbf{v}=\nabla\phi</math>,
 
the energy crossing a surface <math>S</math> in a time period <math>\left[  t,t+T\right]  </math>
 
is
 
<center><math>
 
\rho\int_{t}^{t+T}\int_{S}\phi_{t}\left(  x,y,z,t\right)  \phi_{n}\left(
 
x,y,z,t\right)  d\sigma dt.
 
</math></center>
 
For time-harmonic waves, i.e., <math>\phi=Re \left[  \phi\exp\left(
 
i\omega t\right)  \right]  <math>, and </math>T=2\pi/\omega</math>, the above equation becomes
 
<center><math>
 
-i2\rho\pi\int_{S}\left[  \phi^{}\phi_{n}-\phi\phi_{n}^{}\right]
 
d\sigma=2\rho\piIm \int_{S}\phi^{}\phi_{n}d\sigma
 
. (3-52)
 
</math></center>
 
We notice that if <math>S</math> is the closed surface of domain <math>\mathcal{V}</math>, then from Green's theorem, we have
 
<center><math>
 
\int_{S}\left[  \phi^{}\phi_{n}-\phi\phi_{n}^{}\right]  d\sigma
 
=\int_{\mathcal{V}}\left[  \phi^{}\nabla^{2}\phi-\phi\nabla^{2}\phi^{
 
}\right]  dV=0, (3-82)
 
</math></center>
 
which states that there is no net energy propagation to infinity, i.e., the
 
law of energy conservation. Hence, only the travelling mode <math>H_{0}^{\left(
 
1\right)  }\left(  q_{=T= }r\right)  <math> for the point load and </math>\exp\left(
 
iq_{=T= }r\right)  </math> for the line load carry the energy to infinity since
 
the rest of the solution decays exponentially. Note that the decaying rate of
 
the Hankel function of a real variable is <math>1/\sqrt{r}</math> and the integral on the cylindrical surface gives <math>d\sigma=rdrd\theta dz</math>, thus we have non-zero energy flux for the travelling mode of the point load response.
 
 
 
The evanescent modes in Eqn.~((3-19)) and ((3-18)) are always real since <math>\mathrm{i}R_{n}</math> in Eqn.~((3-18)) and Eqn.~((3-19)) is real for each <math>n=1,2,\cdots</math>. Hence, the only imaginary term in the responses due to point or line load that give non-zero in Eqn.~((3-52)) is the coefficient of the travelling wave. This corresponds to the travelling waves being the only modes that carry energy away from the load. The damped-travelling and evanescent modes contribute motion that is in phase with the forcing, while the travelling mode has a component in quadrature to the forcing.
 
 
 
The velocity potential in the water can be found using Eqn.~((3-5)) to
 
Eqn.~((3-6)), giving
 
<center><math>
 
\phi_{=P= }\left(  r,z\right)  =-\frac{\omega}{2}\sum_{q\in K^{
 
}}\frac{R\left(  q\right)  }{\sinh qH}H_{0}^{\left(  1\right)  }\left(
 
qr\right)  \cosh q\left(  z+H\right)
 
</math></center>
 
for the point load and
 
<center><math>
 
\phi_{=L= }\left(  r,z\right)  =-\omega\sum_{q\in K^{}}
 
\frac{R\left(  q\right)  }{q\sinh qH}\exp\left(  \mathrm{i}qr\right)
 
\cosh q\left(  z+H\right)  .
 
</math></center>
 
for the line load.
 
 
 
Note that the zeros of <math>d\left(  \gamma\right)  \gamma\tanh\gamma H</math>\ are the
 
same as those of the dispersion equation. Thus the singularities of <math>\hat
 
{\phi}<math> and </math>\hat{w}<math> are the same. The term </math>\sinh qH</math> is close to zero for
 
most imaginary roots so these expressions are not directly suitable for
 
numerical computation. The substitution following Eqn.~((3-16)) may be
 
used to give computationally stable expressions.
 
 
 
Each term of <math>w_{=P= }</math> and <math>w_{=L= }</math> in Eqn.~((3-18)) and
 
Eqn.~((3-19)) is a natural mode of a floating ice sheet whose wavenumber
 
is a root of the dispersion equation. We showed that the mode of <math>q_{=T=
 
}</math>, the travelling mode carries the wave energy outwards, and the modes of
 
<math>q_{=D= }</math> and <math>\mathrm{i}q_{n}</math>, <math>n=1,2,...</math>, the damped
 
travelling and evanescent modes, are exponentially decaying.
 
 
 
====Summary of the analytic structure of <math>\hat{w====\left(
 
\gamma\right)  </math>}
 
 
 
Here we summarize and clarify the analytic properties of the complex valued
 
function <math>\hat{w}\left(  \gamma\right)  </math> that have been discussed so far.
 
 
 
It may be worth reminding ourselves that the Mittag-Leffler expansion of <math>\hat{w}</math> given by Eqn.~((ap-20)) has an extra term <math>\hat{w}\left( 0\right)  =0</math> and is unique. Hence, even though a function such as <math>\hat{w}\left(  \gamma\right)  +1</math> has the same singularities and residues, the resulting series expansion will be different. This uniqueness of the expansion guaranties the uniqueness of the solution derived using the inverse Fourier transform of the series expansion form.
 
 
 
We consider the question of whether we can reconstruct the Fourier transform <math>\hat{w}</math> from the modified formula of residues given by Eqn.~((3-17)) and the positions of the roots given by the dispersion equation <math>d\left( \gamma\right)  =0</math>. In addition to the residues and poles, we require that <math>\hat{w}\left(  0\right)  =0</math>, to uniquely reconstruct <math>\hat{w}</math>. The answer to the question is yes and no depending on how the formulae are used. It follows that the function <math>\hat{w}\left(  \gamma\right)  </math> for <math>\gamma\in C</math> can be reconstructed from <math>q\in K^{}</math> and <math>R\left(  q\right)  </math> (either Eqn.~((3-16)) or Eqn.~((3-17))) as
 
<center><math>
 
\hat{w}\left(  \gamma\right)  =\sum_{q\in K^{}}\frac{2qR\left(
 
q\right)  }{\gamma^{2}-q^{2}}.
 
</math></center>
 
However, we find the following
 
<center><math>
 
\hat{w}\left(  \gamma\right)  \neq R\left(  \gamma\right)  \sum_{q\in
 
K^{}}\frac{2q}{\gamma^{2}-q^{2}}. (3-80)
 
</math></center>
 
We note that near <math>\gamma=q\in K^{}</math> the left and the
 
right hand sides share the same analytic property, i.e.,
 
<center><math>
 
\int_{C_{\varepsilon}\left(  q\right)  }\hat{w}\left(  \gamma\right)
 
d\gamma=\int_{C_{\varepsilon}\left(  q\right)  }R\left(  \gamma\right)
 
\sum_{q\in K^{}}\frac{2q}{\gamma^{2}-q^{2}}\,d\gamma
 
</math></center>
 
where <math>C_{\varepsilon}\left(  q\right)  </math> is a circular contour of small radius <math>\varepsilon</math> around a pole <math>q</math>. The reason for ((3-80)) is because the function <math>R\left(  \gamma\right)  </math> defined by Eqn.~((3-16)) and Eqn.~((3-17)) introduce zeros and singularities of their own, which change the analytic properties of the right hand side of Eqn.~((3-80)).
 
 
 
It may seem trivial that the formulae of the residues Eqn.~((3-16)) and Eqn.~((3-17)) are valid only at the pole <math>\gamma=q</math>, nevertheless we have confirmed Eqn.~((3-80)).
 
 
 
===Deep water solution===
 
 
 
It has been mentioned that the imaginary roots lie in the interval
 
<math>\mathrm{i}q_{n}\in\left(  \mathrm{i}\left(  n-1/2\right)  \pi/H,\mathrm{i}n\pi/H\right)  </math> as seen in Fig.~((fig:3-12)). Furthermore the <math>\mathrm{i}q_{n}</math> become equally spaced <math>n\pi/H</math> as the depth <math>H</math> becomes large. Therefore the summation over the evanescent modes is taken at the equally spaced points which will become closely placed as <math>H</math> tends to infinity. This inspires us to find an integral expression of the summation when <math>H=\infty</math>.
 
 
 
We show here that the infinite summations over the evanescent mode of the solutions take particularly simple forms when the water is very deep, i.e., in the limit <math>H\rightarrow\infty</math>. In the limit <math>\left\{  \mathrm{i}q_{n}\right\}  </math> form a continuum with equal density over the imaginary axis. If the residue <math>R\left(  \mathrm{i}q_{n}\right)  </math> decreases proportional to <math>1/H</math>, then the infinite sum over these roots may then be calculated as an integral over the positive imaginary semi-axis. It will be shown that for both types of forcing, the solutions are then given by a sum over special functions at wavenumbers given by the roots of a fifth order polynomial. Furthermore, in the deep-water case <math>H\approx\infty</math>, the solution will be shown to be qualitatively unaffected by setting <math>m=0</math> for typical values of <math>m</math> for sea ice, i.e., smaller than <math>0.1</math>. The governing non-dimensional equations then contain no coefficients depending on the physical properties of the ice or water.
 
 
 
The residue expressed by Eqn.~((3-16)) oscillates more rapidly as the water depth <math>H</math> tends to infinity due to the tangent function, thus the formula becomes unsuitable for numerical computation of the residue. On the other hand, the formula ((3-17)) is smooth and rapidly decaying function as <math>H</math> becomes large. The residue at a root of the dispersion equation in Eqn.~((3-17)) tends to
 
<center><math>\begin{matrix}
 
R\left(  q\right)  & =\frac{1}{H}\frac{\omega^{2}q}{\frac{1}{H}\left\{
 
\omega^{2}\left(  5q^{4}+u\right)  \right\}  +\left(  q^{5}+uq\right)
 
^{2}-\omega^{4}} (eqn:resdeep)\\
 
& \rightarrow\frac{1}{H}Q\left(  q\right)
 
\end{matrix}</math></center>
 
as <math>H\rightarrow\infty</math>, where
 
<center><math>
 
Q\left(  q\right)  =\frac{\omega^{2}q}{\left(  q^{5}+uq\right)  ^{2}
 
-\omega^{4}}. (eqn:Q)
 
 
</math></center>
 
</math></center>
The sum over the imaginary roots in the response to point loading in
+
for <math>Im q>0</math>, <math>r>0</math> can be used. This gives us
Eqn.~((3-18)) is then
 
<center><math>\begin{matrix}
 
\frac{1}{\pi}\sum_{n=1}^{\infty}\mathrm{i}q_{n}R_{n}K_{0}\left(
 
q_{n}r\right)  & \rightarrow\frac{1}{\pi^{2}}\sum_{n=1}^{\infty}
 
\frac{\mathrm{i}n\pi}{H}Q\left(  \frac{\mathrm{i}n\pi}
 
{H}\right)  K_{0}\left(  \frac{n\pi}{H}r\right)  \frac{\pi}{H}\\
 
& \rightarrow\frac{1}{\pi^{2}}\int_{0}^{\infty}\mathrm{i}\xi Q\left(
 
\mathrm{i}\xi\right)  K_{0}\left(  \xi r\right)  dq. (3-62)
 
\end{matrix}</math></center>
 
Similarly, the sum over imaginary roots for line load in Eqn.~((3-19)) is
 
 
<center><math>
 
<center><math>
\mathrm{i}\sum_{n=1}^{\infty}R_{n}\exp\left(  -q_{n}r\right)
+
w\left(  r\right) =\mathrm{i}\sum_{n=-2}^{\infty}
\rightarrow\frac{1}{\pi}\int_{0}^{\infty}\mathrm{i}Q\left(
+
R\left( q\right)  \exp\left(  \mathrm{i}qr\right).
\mathrm{i}\xi\right)  \exp\left(  -\xi r\right) \,d\xi. (3-63)
 
 
</math></center>
 
</math></center>
  
The integrals in Eqn.~((3-62)) and Eqn.~((3-63)) can be evaluated by
+
==Solution for the Velocity Potential==
first writing <math>Q</math> as a sum over simple poles, as we did for the inverse
 
Fourier transform. We notice that <math>Q\left(  q\right)  =\left(  v\left(
 
q\right)  -v\left(  -q\right)  \right)  /2</math> where
 
<center><math>
 
v\left(  \mathrm{i}\xi\right)  =\frac{q}{q^{5}+uq-\omega^{2}
 
} (3-26)
 
</math></center>
 
Then <math>\omega\neq0</math>, <math>v\left(  q\right)  </math> and <math>v\left(  -q\right)  </math> have no
 
poles in common and hence poles of <math>v</math> are also poles of <math>Q</math>. It follows that
 
all poles of <math>Q</math> are plus and minus the poles of <math>v</math>. Since <math>v</math> and <math>w</math> (in
 
the limit <math>H\rightarrow\infty</math>) coincide for arguments with positive real
 
part, the poles of <math>\hat{w}</math> with positive real part, i.e. <math>q_{=T=
 
},q_{=D= }<math>, and </math>q_{=D= }^{}<math>, are also roots of </math>v</math> and hence
 
<math>Q</math>. There are two further roots of <math>v</math> with negative real part, which we
 
denote <math>q_{=E= }</math> and <math>q_{=E= }^{}</math>, with <math>q_{=E= }</math> chosen to
 
have positive imaginary part. It can be shown that <math>Im \left(
 
q_{=E= }\right)  >0<math> for any finite mass density </math>m</math>, and hence
 
<math>q_{=E= } </math> is never real.
 
 
 
\begin{figure}[tbh]
 
\begin{center}
 
\includegraphics[height=5.2455cm,width=7.0534cm]{deeprt.eps}
 
\caption{Relative positions of the poles for the deep-water response. <math></math>,<math>\triangle</math> and <math>\diamondsuit</math> indicate <math>q_{=T= }</math>, <math>q_{=D= }</math> (<math>q_{=D= }^{}</math>) and <math>q_{=E= }</math> (<math>q_{=E= }^{}</math>).}
 
(fig:3-25)
 
\end{center}
 
\end{figure}
 
  
Let <math>K_{v}=\left\{  q_{=T= },q_{=D= },q_{=D= }^{},q_{=E=  },q_{=E= }^{}\right\}  ,</math> shown in Fig.~((fig:3-25)), denote the set of poles of <math>v</math>. The residue of <math>v</math> at a pole <math>q\in K_{v}</math> is
+
The velocity potential in the water can be found in two-dimensions
 
<center><math>
 
<center><math>
R_{v}\left(  q\right)  =\frac{q^{2}}{5\omega^{2}-4uq}. (eqn:resv)
+
\phi\left(  r,z\right)  =-\frac{\omega}{2}\sum_{n=-2}^{\infty}
 +
\frac{R\left(  q_n\right)  }{\sinh q_nH}H_{0}^{\left(  1\right)  }\left(
 +
q_nr\right)  \cosh q_n\left( z+H\right)
 
</math></center>
 
</math></center>
The residues at <math>q_{=T= }</math>, <math>q_{=D= }</math> and <math>q_{=D= }
+
and for three-dimensions
^{*}</math>, are the same as defined in section 3.2.4. We respectively denote these
 
<math>R_{=T= }</math>, <math>R_{=D= }</math> and <math>R_{=D= }^{}</math>, as before,
 
respectively. Write <math>R_{=E= }=R_{v}\left(  q_{=E= }\right)  </math> and
 
hence <math>R_{=E= }^{}=R_{v}\left(  q_{=E= }^{}\right)  </math>. The
 
residues of <math>Q</math> are <math>1/2</math> the residue of <math>v</math>.
 
 
 
The integral in Eqn.~((3-62)) may be evaluated using the fractional function expansion
 
 
<center><math>
 
<center><math>
\mathrm{i}\xi Q\left(  \mathrm{i}\xi\right)  =-\sum_{q\in K_{v} }\frac{q^{2}R_{v}\left(  q\right)  }{\xi^{2}+q^{2}} (eqn:ikQexpn)
+
\phi\left(  r,z\right)  =-\omega\sum_{n=-2}^{\infty}
 +
\frac{R\left(  q_n\right)  }{q_n\sinh q_nH}\exp\left(  \mathrm{i}q_nr\right)
 +
\cosh q_n\left( z+H\right).
 
</math></center>
 
</math></center>
and the integral (Abramowitz and Stegun [[Stegun]] formula 11.4.47 with
 
<math>\nu=0</math>, <math>t=k</math>, <math>r=a</math>),
 
<center><math>
 
\int_{0}^{\infty}\frac{K_{0}\left(  \xi r\right)  }{\left(  \xi^{2}
 
+q^{2}\right)  }d\xi=\frac{\pi^{2}}{4q}\left[  \mathbf{H}_{0}\left(
 
qr\right)  -Y_{0}\left(  qr\right)  \right]  ,
 
</math></center>
 
which holds for <math>Re q>0</math> since we always take <math>r>0</math>. Here <math>\mathbf{H}_{0}</math> is a Struve function of zero order and its power expansion for numerical computation is shown in appendix A. For notational brevity, we use the function <math>h\left(  rz\right)  =\mathbf{H}_{0}\left(  rz\right) -Y_{0}\left(  rz\right)  </math>. Performing the integral in Eqn.~((3-62)) and combining conjugate pair poles, Eqn.~((3-18)) for the response to point load in the deep water limit can be written in terms of poles of <math>v</math> as,
 
<center><math>\begin{matrix}
 
w_{=P= }\left(  r\right)  & =\frac{\mathrm{i}}{2}q_{=T=
 
}R_{=T= }H_{0}^{\left(  1\right)  }\left(  q_{=T= }r\right)
 
-Im \left[  q_{=D= }R_{=D= }H_{0}^{\left(  1\right)
 
}\left(  q_{=D= }r\right)  \right] \\
 
& -\frac{q_{=T= }R_{=T= }}{4}h\left(  q_{=T= }r\right)  -\frac
 
{1}{2}Re \left(  q_{=D= }R_{=D= }h\left(  q_{=D=
 
}r\right)  \right)  +\frac{1}{2}Re \left(  q_{=E=
 
}R_{=E= }h\left(  -q_{=E= }r\right)  \right)  . (3-60)
 
\end{matrix}</math></center>
 
 
The pole <math>-q_{=E= }</math> of <math>Q</math> has been used since <math>Re \left( -q_{=E= }\right)  >0</math>. Note that only the first, travelling wave term is pure imaginary number, which corresponds to that mode being the only one that propagates energy away from the point of forcing. The remaining terms give displacements that are in phase with the forcing.
 
 
The integral in Eqn.~((3-63)) for line load may be found using the expansion
 
<center><math>
 
\mathrm{i}Q\left(  \mathrm{i}\xi\right)  =\sum_{q\in K_{v}}
 
\frac{\xi R_{v}\left(  q\right)  }{\xi^{2}+q^{2}}
 
</math></center>
 
and the integral (Abramowitz and Stegun [[Stegun]] formula 5.2.13 with
 
<math>t=k/q</math> and hence <math>z=qr</math> and formula 5.2.7)
 
<center><math>
 
\int_{0}^{\infty}\frac{\xi\exp\left(  -\xi r\right)  }{\xi^{2}+q^{2}}
 
d\xi=-Ci \left(  qr\right)  \cos\left(  qr\right)
 
-=si= \left(  qr\right)  \sin\left(  qr\right)
 
</math></center>
 
holding for <math>Re \left(  q\right)  >0</math> since <math>r</math> is positive
 
real. Here <math>Ci </math> and <math>si </math> are cosine integral
 
and sine integral functions (see appendix A), respectively. As with the
 
point-forcing case, the notation and computation are simplified by defining
 
the function <math>g\left(  qr\right)  =-Ci \left(  qr\right)
 
\cos\left(  qr\right)  -si \left(  qr\right)  \sin\left(
 
qr\right)  </math>. Then, combining conjugate-pair poles, Eqn.~((3-19)) for the
 
response to line forcing in the deep water limit can be written as,
 
<center><math>\begin{matrix}
 
w_{=L= }\left(  r\right)  & =\mathrm{i}R_{=T= }\exp\left(
 
\mathrm{i}q_{=T= }r\right)  -2Im \left[  R_{=D=
 
}\exp\left(  \mathrm{i}q_{=D= }r\right)  \right] \\
 
& +\frac{R_{=T= }}{\pi}g\left(  q_{=T= }r\right)  +\frac{2}{\pi
 
}Re \left[  R_{=D= }g\left(  q_{=D= }r\right)  \right]
 
+\frac{2}{\pi}Re \left[  R_{=E= }g\left(  -k_{=E=
 
}r\right)  \right]  . (3-61)
 
\end{matrix}</math></center>
 
where, again, the pole <math>-q_{=E= }</math> has been used.
 
 
We emphasize here that the denominator of the function <math>v</math> in Eqn.~((3-26)) is not meant to be the dispersion function of the deep water problem, but it is only the partial expression of the dispersion function. The full deep-water dispersion equation may be written as
 
<center><math>
 
\gamma^{4}-m\omega^{2}+1-\frac{\omega^{2}}{\gammasgn \left(
 
Re \gamma\right)  }=0 (3-28)
 
</math></center>
 
where <math>sgn \left(  Re \gamma\right)  </math> is
 
sign-function of the real part of <math>\gamma</math>, which has continuous singularity
 
on the imaginary axis. We have obtained the deep-water solution without
 
dealing with this continuous singularity. The function <math>\gamma
 
sgn \left(  Re \gamma\right)  </math> is often defined
 
using the limit of <math>\sqrt{\gamma^{2}\pm i\varepsilon^{2}}</math> as <math>\varepsilon
 
\rightarrow0</math>, which is defined on a two-sheeted Riemann surface and denoted
 
by <math>\left|  \gamma\right|  _{\pm}</math> because it is equal to <math>\left|
 
\gamma\right|  </math> on the real axis.
 
 
===Computation of the solutions===
 
 
====Static load====
 
 
We can find the displacement <math>w_{=P= }</math> and <math>w_{=L= }</math>\ for static
 
loading by setting <math>\omega=0</math> in Eqn.~((3-9)), to give
 
<center><math>
 
d\left(  \gamma\right)  =\gamma^{4}+1 (3-35)
 
</math></center>
 
and following the same procedure shown in section 3.2. The Fourier transform
 
of the displacement is
 
<center><math>
 
\hat{w}\left(  \gamma\right)  =\frac{1}{\gamma^{4}+1}
 
</math></center>
 
which has four poles, <math>\pm e^{\mathrm{i}\pi/4}</math>, <math>\pm
 
e^{\mathrm{i}3\pi/4}<math>, the set of poles in the upper half plane being </math>K^{}=\left\{  e^{\mathrm{i}\pi/4},e^{\mathrm{i} 3\pi/4}\right\}  <math>. The residue of </math>\hat{w}\left(  \gamma\right)  <math> at a pole </math>q</math> is
 
<center><math>
 
R\left(  \gamma\right)  =\left.  \frac{1}{\left(  \gamma^{4}+1\right)
 
^{\prime}}\right|  _{\gamma=q}=\frac{1}{4q^{3}}.
 
</math></center>
 
Hence, the solution is, from Eqn.~((3-14))
 
<center><math>\begin{matrix}
 
w_{=P= }\left(  r\right)  & =\frac{1}{4\pi}\left\{  e^{-\mathrm{i} \pi/2}K_{0}\left(  -\mathrm{i}e^{\mathrm{i}\pi/4}r\right)
 
+e^{-\mathrm{i}3\pi/2}K_{0}\left(  -\mathrm{i}e^{\mathrm{i}3\pi/4}r\right)  \right\} \\
 
& =\frac{\mathrm{i}}{4\pi}\left\{  -K_{0}\left(  e^{-\mathrm{i} \pi/4}r\right)  +K_{0}\left(  e^{\mathrm{i}\pi/4}r\right)  \right\} \\
 
& =-\frac{kei \left(  r\right)  }{2\pi} (3-32)
 
\end{matrix}</math></center>
 
where <math>kei \left(  \zeta\right)  </math> is the Kelvin function (of zero order) and we have used the identity (Abramowitz and Stegun [[Stegun]] formulae 9.9.2 and 9.6.32)
 
<center><math>
 
\mathrm{i}2kei \left(  \zeta\right)  =K_{0}\left(
 
e^{\mathrm{i}\pi/4}\zeta\right)  -K_{0}\left(  e^{-\mathrm{i}
 
\pi/4}\zeta\right)  .
 
</math></center>
 
The solution given by Eqn.~((3-32)) is the same solution derived by
 
Wyman [[Wyman50]].
 
Using Eqn.~((3-15)), we can derive the static line-loading solution
 
<center><math>
 
w_{=L= }\left(  r\right)  =\frac{1}{2}\exp\left(  \frac{-r}{\sqrt{2}
 
}\right)
 
</math></center>
 
 
====Deflection at the location of forcing====
 
 
It appears that expression ((3-18)) has a singularity at the origin since
 
the Hankel function and the modified Bessel function behave as <math>\log</math>-like
 
function near the origin. Indeed, Strathdee, Robinson and Haines [[Strathdee91]]
 
erroneously stated that the solution for a floating thin plate had a
 
singularity at the point of forcing. However, since the plate equation is
 
fourth order and Laplace's equation is second order, we would expect the
 
solution to be smooth everywhere, including at <math>r=0</math>. We show here that the
 
summation is indeed bounded everywhere and derive an expression for the
 
displacement at the point of forcing that is convenient to compute.
 
 
The modified Bessel function <math>K_{0}\left(  \zeta\right)  </math> has the polynomial
 
form
 
<center><math>
 
K_{0}\left(  \zeta\right)  =-\log\left(  \frac{\zeta}{2}\right)  I_{0}\left(
 
\zeta\right)  +\sum_{l=0}^{\infty}\frac{\left(  \zeta/2\right)  ^{2l}}{\left(
 
l!\right)  ^{2}}\psi\left(  l+1\right)  ,
 
</math></center>
 
where <math>I_{0}</math> and <math>\psi</math> are the modified Bessel function and the Psi
 
function, respectively (Abramowitz and Stegun [[Stegun]] formulae 9.6.12
 
and 9.6.13, appendix A). For small <math>\left|  \zeta\right|  </math>, <math>K_{0}\left(
 
\zeta\right)  \approx-\log\zeta+c<math> where </math>c=\log2+\psi\left(  1\right)  </math>
 
since <math>I_{0}\left(  0\right)  =1</math>. Hence, as <math>r\rightarrow0</math>, the infinite
 
series of <math>w_{=P= }</math> becomes a series of log-function of <math>r</math>,
 
<center><math>\begin{matrix}
 
w_{=P= }\left(  r\right)  & \rightarrow\frac{1}{\pi}\sum_{q\in
 
K^{}}qR\left(  q\right)  \left(  -\log\left(  -\mathrm{i}
 
qr\right)  +c\right) \\
 
& =-\frac{1}{\pi}\sum_{q\in K^{}}qR\left(  q\right)  \log\left(
 
-\mathrm{i}q\right)  +\frac{c-\log r}{\pi}\sum_{q\in K^{}
 
}qR\left(  q\right)  . (3-42)
 
\end{matrix}</math></center>
 
We notice that the second term becomes singular as <math>r\rightarrow0</math>.
 
 
Consider now a contour integration of the function <math>\hat{w}\left(
 
\gamma\right)  \gamma</math>, anti-clockwise along the contour shown in
 
Fig.~((fig:3-15)). The arc of radius <math>A</math> is chosen to avoid the poles on
 
the imaginary axis, and the arcs around the poles on the real-axis are taken
 
to have small radius.
 
 
\begin{figure}[tbh]
 
\begin{center}
 
\includegraphics[height=4.0857cm,width=6.4647cm]{contour.eps}
 
\caption{Contour used for integration with approximate pole positions shown as crosses.}
 
(fig:3-15)
 
\end{center}
 
\end{figure}
 
 
Since <math>\hat{w}\left(  \gamma\right)  \gamma</math> is an odd function, the integral
 
over the real axis, including the two small arcs, is zero. Further, as
 
<math>A\rightarrow\infty</math>, <math>\hat{w}\left(  \gamma\right)  \gamma</math> tends to zero
 
faster than <math>A^{-2}</math> on the semi-circle of radius <math>A</math>, and the integral over
 
the semi-circle tends to zero as <math>A\rightarrow\infty</math>. Hence the integral over the whole contour tends to zero as <math>A\rightarrow\infty</math>. Since this limit
 
equals a constant multiplied by the sum of the residues of <math>\hat{w}\left(
 
\gamma\right)  \gamma</math> at the poles enclosed within the contour, the sum of
 
residues of <math>\hat{w}\left(  \gamma\right)  \gamma</math> at poles in the upper half
 
plane is zero, i.e.,
 
<center><math>
 
0=\int_{C}\hat{w}\left(  \gamma\right)  \gamma d\gamma=2\pi\mathrm{i}
 
\sum_{q\in K^{}}qR\left(  q\right)  .
 
</math></center>
 
We immediately see that the term multiplied by <math>\left(  c-\log r\right)  /\pi</math> in Eqn.~((3-42)), which is the singular part, is zero. Thus at <math>r=0</math> the complex displacement takes the finite value
 
<center><math>
 
w_{=P= }\left(  0\right)  =-\frac{1}{\pi}\sum_{q\in K^{}
 
}qR\left(  q\right)  \log\left(  -\mathrm{i}q\right)  =-\frac{1}{\pi
 
}\sum_{q\in K^{}}qR\left(  q\right)  \log\left(  q\right)
 
. (3-45)
 
</math></center>
 
Since <math>qR\left(  q\right)  </math> decreases as <math>q_{n}^{-8}\propto n^{-8}</math> for the
 
evanescent modes, relatively few terms are required to evaluate this sum
 
accurately. Again, the computation of the solution requires fewer modes
 
when\ the depth <math>H</math> is small because of the same reason mentioned earlier.
 
 
Fig.~((fig:3-26)) shows the displacement at the point of forcing as a
 
function of frequency and for water depths <math>H=20\pi</math>, <math>2\pi</math>, and <math>0.2\pi</math>.
 
The depths are taken as multiples of <math>2\pi</math> since that is the wavelength of
 
the travelling wave with unit non-dimensional wavenumber.
 
 
\begin{figure}[tbh]
 
\begin{center}
 
\includegraphics[height=9.8079cm,width=11.3368cm]{pload.eps}
 
\caption{Top to bottom shows the magnitude and argument part of the complex
 
displacement at the point of forcing as a function of non-dimensional
 
frequency. The non-dimensional water depths <math>H=2\pi\times10</math>, 2<math>\pi\times1</math>,
 
2<math>\pi\times0.1</math>, are shown. In all cases we take <math>m=0</math>.}
 
(fig:3-26)
 
\end{center}
 
\end{figure}
 
 
Water depths greater than <math>20\pi</math> give visually identical deflections, so the
 
curve for <math>H=20\pi</math> may be taken as the deep-water solution. That is the
 
solution for <math>H=2\pi</math> is nearly identical to the deep-water solution at all
 
frequencies, so that <math>H=2\pi</math> may be considered ''deep''  for point load.
 
 
We notice by spliting Eqn.~((3-45)) into the real and imaginary parts,
 
<center><math>\begin{matrix}
 
w_{=P= }\left(  0\right)  & =-\frac{1}{\pi}\left(  q_{=T=
 
}R_{=T= }\log\left(  q_{=T= }\right)  +2Re \left[
 
q_{=D= }R_{=D= }\log\left(  -\mathrm{i}q_{=D= }\right)
 
\right]  +\sum_{n=1}^{\infty}\mathrm{i}q_{n}R_{n}\log\left(
 
q_{n}\right)  \right) \\
 
& +\mathrm{i}\frac{q_{=T= }R_{=T= }}{2}.
 
\end{matrix}</math></center>
 
The imaginary part is just the coefficient of the travelling mode. This
 
corresponds to the travelling mode being the only mode that carries energy
 
away from the point of forcing.
 
 
The displacement at the point of forcing when the ice is floating on deep water may be found by evaluating the finite-depth solution in Eqn.~((3-45)) in the limit <math>H\rightarrow\infty</math> as described in the previous section. The sum over evanescent modes is
 
<center><math>
 
-\frac{1}{\pi}\sum_{n=1}^{\infty}\mathrm{i}q_{n}R\left(
 
\mathrm{i}q_{n}\right)  \log\left(  q_{n}\right)  \rightarrow-\frac
 
{1}{\pi^{2}}\int_{0}^{\infty}\mathrm{i}qQ\left(  \mathrm{i}
 
q\right)  \log\left(  q\right)  dq.
 
</math></center>
 
Using the expansion in Eqn.~((eqn:ikQexpn)) and the integral identity
 
(Erdelyi [[bateman2]] section 14.2 (24) with <math>a=-\mathrm{i}q</math>, and
 
<math>y=\mathrm{i}q</math>)
 
<center><math>
 
\int_{0}^{\infty}\frac{\log x}{x^{2}+q^{2}}dx=\frac{\pi}{2q}\log q
 
</math></center>
 
for <math>Re q>0</math>,\ we find that
 
<center><math>
 
w_{=P= }\left(  0\right)  =\mathrm{i}\frac{q_{=T= }R_{=T=
 
}}{2}-\frac{q_{=T= }R_{=T= }}{2\pi}\log\left(  q_{=T= }\right)
 
-\frac{1}{\pi}Re \left(  q_{=D= }R_{=D= }\log\left(
 
-q_{=D= }\right)  \right)  -\frac{1}{\pi}Re \left(
 
q_{=E= }R_{=E= }\log\left(  -q_{=E= }\right)  \right)
 
. (3-41)
 
</math></center>
 
This gives a very simple expression for finding the displacement at the point of forcing when the water is deep. Hence, Eqn.~((3-41)) gives an efficient route to the result computed numerically by Nevel [[Nevel70]].
 
 
Because the exponentials in Eqn.~((3-19)) are continuous everywhere, the displacement at the line of forcing may found by setting <math>r=0 </math> directly to give
 
<center><math>
 
w_{=L= }\left(  0\right)  =\sum_{q\in K^{}}\mathrm{i}
 
R\left(  q\right)  . (3-46)
 
</math></center>
 
Using this expression, the absolute value, and argument parts of the complex
 
deflection, plotted as a function of frequency, are shown in
 
Fig.~((fig:3-16)) for water depths of <math>H=20\pi</math>, <math>2\pi</math>, and <math>0.2\pi</math>.
 
 
\begin{figure}[tbh]
 
\begin{center}
 
\includegraphics[height=9.7069cm,width=11.2621cm]{lload.eps}
 
\caption{Top to bottom shows the magnitude, argument part of the complex displacement on the line of forcing as a function of non-dimensional frequency for the non-dimensional water depths <math>H=2\pi\times10</math>, <math>2\pi\times1</math>, <math>2\pi\times0.1</math>. We take <math>m=0</math>.}
 
(fig:3-16)
 
\end{center}
 
\end{figure}
 
 
We see little difference in Fig.~((fig:3-16)) between the graphs for
 
<math>H=2\pi</math> and <math>20\pi</math>, thus <math>H</math> greater than <math>2\pi</math> may be considered as
 
''deep''  for line forcing. However, at the frequency lower than <math>1</math> the
 
graphs show difference between <math>H=2\pi</math> and <math>20\pi</math>, hence <math>H=2\pi</math> may not be considered as ''deep''  at low frequency for line forcing.
 
 
Because the cosine integral has a log-like singularity at the origin,
 
Eqn.~((3-61)) is not directly suitable for computing the displacement on
 
the line of forcing in the deep-water limit. However, expanding
 
<math>Ci \left(  qr\right)  </math> in terms of the log plus power series
 
and using the identity <math>\sum_{q\in K_{v}}R_{v}\left(  q\right)  =0</math> allows us
 
to write
 
<center><math>
 
w_{=L= }\left(  0\right)  =\mathrm{i}R_{=T= }-\frac
 
{R_{=T= }}{\pi}=log= \left(  w\right)  -\frac{2}{\pi}Re
 
\left[  R_{=D= }=log= \left(  -q_{=D= }\right)  \right]  -\frac
 
{2}{\pi}Re \left[  R_{=E= }=log= \left(  -q_{=E=
 
}\right)  \right].  (3-64)
 
</math></center>
 
The above formula is easily evaluated to give the displacement on the line of forcing when the water is deep.
 
 
====Derivatives of the deflection====
 
 
As mentioned, the modified Bessel function <math>K_{0}</math> and Hankel function
 
<math>H_{0}^{\left(  1\right)  }</math> have a log-like singularity at <math>r=0</math>. Hence
 
depending on the software package used for numerical computation, the
 
solutions may not be stable near the origin. In this section we show
 
alternative formulae for the derivatives of the solution that is stable near
 
<math>r=0</math> using a power series expansion of <math>I_{0}</math> given in Abramowitz and Stegun [[Stegun]] formula 9.6.10 with <math>\nu=0</math>. We write the deflection as the sum of the deflection at <math>r=0</math> and the remainder,
 
<center><math>
 
w_{=P= }\left(  r\right)  =w_{=P= }\left(  0\right)  +\frac{1}{\pi
 
}\sum_{q\in K^{}}qR\left(  q\right)  \left(  \sum_{l=1}^{\infty
 
}\frac{\left(  -q^{2}r^{2}\right)  ^{l}}{4^{l}\left(  l!\right)  ^{2}}\left(
 
\psi\left(  l+1\right)  -\log\left(  \frac{-\mathrm{i}qr}{2}\right)
 
\right)  \right)  (3-65)
 
</math></center>
 
for all <math>r</math>. Note that near the origin the deflection behaves as <math>\left(
 
\log\left(  qr\right)  +c\right)  r^{2}</math> and thus the first derivative of the
 
deflection near the origin is
 
<center><math>
 
\frac{d}{dr}\left(  \left(  \log\left(  qr\right)  +c\right)  r^{2}\right)
 
=r\left(  1+2\left(  \log\left(  qr\right)  +c\right)  \right)  .
 
</math></center>
 
Therefore, we have <math>w_{=P= }^{\prime}\left(  0\right)  =0</math> as expected. It follows that the displacement function obtained by the method above is regular everywhere.
 
 
Measurements of flexural response of ice sheets can be made using the strain
 
gauge which measures the curvature of the upper surface of the ice sheet.
 
Measuring the strain will be examined further in section 3.7. The second
 
derivative with respect to <math>r</math> for each term in Eqn.~((3-65)) with
 
<math>l\geq1</math> has the <math>r</math>-dependence
 
<center><math>
 
\frac{d^{2}}{dr^{2}}\left(  \left(  \log\left(  qr\right)  +c\right)
 
r^{2l}\right)  =r^{2l-2}\left(  \log\left(  qr\right)  \left(  4l^{2}
 
-2l\right)  +4l^{2}c+4l-2lc-1\right)
 
</math></center>
 
which is non-zero at <math>r=0</math> only for the <math>l=1</math> term. Since <math>\sum_{q\in
 
K^{}}R\left(  q\right)  q^{3}=1/2</math>, evaluating the integral of
 
<math>\hat{w}\left(  \gamma\right)  \gamma^{3}</math> on the contour shown in
 
Fig.~((fig:3-15)), the second derivative of the <math>l=1</math> term in
 
Eqn.~((3-65)) can be written
 
<center><math>
 
\frac{d^{2}}{dr^{2}}w_{=P= }\left(  r\right)  \approx\frac{1}{8\pi}\left(
 
2\log\left(  \frac{r}{2}\right)  +3-\psi\left(  2\right)  +\sum_{q\in
 
K^{}}q^{3}R\left(  q\right)  \log\left(  -\mathrm{i}
 
q\right)  \right)
 
</math></center>
 
for small <math>r</math>. So we see that the strain has a singularity at the origin that
 
behaves like <math>\log r</math> and is in phase with the forcing. Using formulae for the
 
derivatives (recursive formulae in Erdelyi [[bateman2]] formula) of the
 
Hankel function, we have the first and second derivatives,
 
<center><math>\begin{matrix}
 
\frac{d}{dr}w_{=P= }\left(  r\right)  & =\sum_{q\in K^{}}
 
q^{2}R\left(  q\right)  H_{1}^{\left(  1\right)  }\left(  qr\right)
 
, (3-71)\\
 
\frac{d^{2}}{dr^{2}}w_{=P= }\left(  r\right)  & =\sum_{q\in
 
K^{}}q^{2}R\left(  q\right)  \left\{  qH_{2}^{\left(  1\right)
 
}\left(  qr\right)  -\frac{1}{r}H_{1}^{\left(  1\right)  }\left(  qr\right)
 
\right\}  . (3-72)
 
\end{matrix}</math></center>
 
 
We may use simpler formulae to compute the deflection using asymptotic
 
formulae of Hankel function of the solutions. It is shown in Abramowitz and
 
Stegun ([[Stegun]] formula 9.2.3) that
 
<center><math>
 
H_{0}^{\left(  1\right)  }\left(  \zeta\right)  \sim\sqrt{\frac{2}{\pi\zeta}
 
}\exp\left\{  \mathrm{i}\left(  \zeta-\pi/4\right)  \right\}
 
</math></center>
 
for large <math>\left|  \zeta\right|  </math>. Because the terms other than the
 
travelling mode decay exponentially, selecting just the term due to the
 
travelling mode in Eqn.~((3-60)) gives the displacement at far field as
 
<center><math>
 
w_{=P= }\left(  r\right)  \sim\frac{\sqrt{q_{=T= }}R_{=T= }}
 
{\sqrt{2\pi}}\frac{\exp\left\{  \mathrm{i}\left(  q_{=T= }
 
r+\pi/4\right)  \right\}  }{\sqrt{r}} (eqn:etapfar)
 
</math></center>
 
for large <math>r</math>.
 
 
\begin{figure} [tbh]
 
\begin{center}
 
\includegraphics[height=7.0973cm,width=10.2297cm]{maxamp.eps}
 
\caption{Curves of the coefficients of <math>w_{=P= }</math> (solid) and <math>w_{=L= }</math> (dash-dot) in the far field given by Eqn. ((eqn:etapfar)) and Eqn. ((eqn:etaplfar)) respectively.}
 
(fig:3-30)
 
\end{center}
 
\end{figure}
 
 
Equivalently, we find from Eqn.~((3-19)) and Eqn.~((3-61)) that the deflection far from a line load is
 
<center><math>
 
w_{=L= }\left(  x\right)  \sim R_{=T= }\exp\left\{  \mathrm{i}
 
\left(  q_{=T= }\left|  x\right|  +\pi/2\right)  \right\}
 
(eqn:etaplfar)
 
</math></center>
 
for large <math>\left|  x\right|  </math>. Fig.~((fig:3-30)) shows the coefficients
 
of Eqn.~((eqn:etapfar)) and Eqn.~((eqn:etaplfar)). The maximum
 
response of point load at the far field is achieved at <math>\omega=0.81</math> and that
 
of line forcing at <math>\omega=0.73</math>.
 
 
\begin{figure}[ptb]
 
\begin{center}
 
\includegraphics[height=12.9689cm,width=11.6992cm]{displ.eps}
 
\caption{Amplitude and phase of the deep-water displacement function <math>w\left( r,\omega\right)  </math> of point load at various non-dimensional frequencies <math>\omega=0.2,1.0 </math> and <math>5.0</math> from the top respectively.}
 
(fig:3-13)
 
\end{center}
 
\end{figure}
 
 
\begin{figure}[ptb]
 
\begin{center}
 
\includegraphics[height=11.9035cm,width=13.0589cm]{dispp.eps}
 
\caption{Amplitude and phase of deep-water displacement function <math>w\left( r,\omega\right)  </math> for the line load problem at various non-dimensional frequencies <math>\omega=0.2,1.0</math> and <math>5.0</math>. The displacement at <math>\omega=0.2</math> for <math>H=2\pi</math> is plotted, since at a higher frequency, the response is indistinguishable between <math>H=2\pi</math> and <math>\infty</math>.}
 
(fig:3-14)
 
\end{center}
 
\end{figure}
 
 
All graphs of the solutions shown in Fig.~((fig:3-13)) and Fig.~((fig:3-14)) are generated using MatLab, and the computer codes are direct implementation of the formulae reported here. The special functions are computed using the built-in functions of MatLab, which are stable enough for small <math>r</math>. The direct formula given by Eqn.~((3-72)) is also stable enough for small <math>r</math> and graphs of the strain are shown in section 3.6. However it is also possible to compute them using the power expansions of the special functions, although that computation is not as fast as the built-in functions.
 
 
Curves in Fig.~((fig:3-13)) show deflections for various non-dimensional radial frequencies <math>\omega=0.2,1.0,5.0</math> for the deep-water case. At low frequency <math>\omega=0.2</math> the deflection is nearly identical to the static solution having the imaginary part of the deflection nearly zero. Curves in Fig.~((fig:3-14)) show comparison of the deflection for a line-loading case between the non-dimensional water depth <math>H=2\pi</math> and <math>\infty</math> at various non-dimensional radial frequencies. The effects of the water depth can be seen only for the low frequency case as expected from the deflection at the origin in Fig.~((fig:3-16)).
 
 
===Scaling of the solution (sec:2)===
 
 
An advantage of scaling or non-dimensionalization of the solution obviously is
 
simplification of the system of equations. However, as shown by Fox
 
[[Fox00]], an appropriate scaling method enables us to find a solution
 
which represents all physical solutions for ranges of physical parameters, so
 
that we only need to obtain the non-dimensional solution and the
 
characteristic length in order to compute all physical solutions. We have not
 
given theoretical justifications for our particular choice of the
 
characteristic length and time, other than the fact that the scaling factors
 
seemed well suited for the plate equation and the dispersion equation. It may
 
seem illogical to say that appropriate scaling factors can only be found after
 
the system of equations is solved and analytical solutions are derived. A
 
properly scaled solution gives us informations\ that are universal to all
 
physical solutions, thus an obvious application may be in scaled model
 
experiments, where, for example, one might interpret scaled tank experiments
 
to the actual size measurements. It is shown by Fox [[Fox00]] that despite
 
the fact that the scaling constants <math>l_{=c= }</math> and <math>t_{=c= }</math> are
 
originally obtained from the problem of static loading of an infinite plate,
 
the same scaling relations hold for an ice sheet of more general shape.
 
 
====Scaled solutions and physical solutions====
 
 
We again consider the solution for the static load shown in subsection 3.5.1.
 
If we solve the original system of differential equations, we find that the
 
dispersion equation is
 
<center><math>
 
d\left(  \bar{\gamma}\right)  =D\bar{\gamma}^{4}+\rho g (3-30)
 
</math></center>
 
and the physical solution is
 
<center><math>
 
\bar{w}_{=P= }\left(  \bar{r}\right)  =-\frac{kei \left(
 
\bar{r}\right)  }{2\pi\sqrt{D\rho g}}=-\frac{kei \left(
 
\bar{r}\right)  }{2\pi\rho gl_{=c= }^{2}} (3-31)
 
</math></center>
 
where <math>\bar{r}</math> is the physical distance from the forcing point. Observing
 
Eqn.~((3-30)) and Eqn.~((3-31)), we find the appropriate
 
non-dimensionalization constant <math>l_{=c= }</math> to be <math>\left(  D/\rho g\right)
 
^{1/4}</math> which of course gave the non-dimensional solution ((3-32)) and
 
then the relationship between the non-dimensional and physical solutions. The
 
conversion between the scaled and physical solutions is found by back-stepping
 
the non-dimensionalization and the Fourier transforms (in one and two
 
dimensional space) shown in section 3.2. Using the notation with the bar for
 
the dimensional variables, we have the following conversion relationship for
 
the unit point load, <math>1</math> Newton
 
<center><math>
 
\bar{w}_{=P= }\left(  \bar{r}\right)  =\frac{1}{\rho gl_{=c= }^{2}
 
}w_{=P= }\left(  r\right)  .
 
</math></center>
 
Note that <math>\bar{w}_{=P= }\left(  \bar{r}\right)  </math> is displacement per
 
Newton.
 
 
Again, considering the scaled and physical solutions of the static
 
line-loading problem, we find the conversion relation. The sum in
 
Eqn.~((3-64)) gives the physical response to a static line force as
 
<center><math>
 
\bar{w}_{=L= }\left(  \left|  \bar{x}\right|  \right)  =\frac{1}{2\rho
 
gl_{=c= }}\exp\left(  \frac{-\left|  \bar{x}\right|  }{\sqrt
 
{2}l_{=c= }}\right)  \cos\left(  \frac{\left|  \bar{x}\right|  }{\sqrt
 
{2}l_{=c= }}-\frac{\pi}{4}\right)  .
 
</math></center>
 
Therefore, for the line load the physical displacement for the line forcing,
 
we have
 
<center><math>
 
\bar{w}_{=L= }\left(  \left|  \bar{x}\right|  \right)  =\frac{1}{\rho
 
gl_{=c= }}w_{=L= }\left(  \left|  x\right|  \right)  .
 
</math></center>
 
 
The physical dimensions are completely absent in the static-load dispersion
 
function ((3-35)). Hence the roots of the dispersion equation are
 
independent of any change in physical parameters, such as ice thickness, water
 
depth, or mass density. In other words, once the formula ((3-32)) is
 
obtained, the characteristic length <math>l_{=c= }</math>\ is the only parameter that
 
we need to define the characteristics of an ice sheet. The expression in
 
Eqn.~((3-32)) is the same as the solution given by Wyman [[Wyman50]].
 
 
One may hope that the same characteristic length <math>l_{=c= }</math>\ and
 
normalizing scheme could be applied to the dynamic ice sheet problem.
 
Fortunately, we find that the same characteristic length <math>l_{=c= }</math>\ and
 
time <math>t_{=c= }</math> may be used to scale the system of equations for the range
 
of ice thickness <math>h</math> and forcing frequency <math>\omega</math>\ that is relevant to the
 
study of sea ice. The validity of the non-dimensionalization of
 
Eqn.~((2-74)) is based primarily on the dispersion equation ((3-28))
 
for deep water. Fig.~((fig:3-5)) shows the positions of the wave-numbers
 
<math>q_{=T= }</math>, <math>q_{=D= }</math> and <math>q_{=E= }</math> of the normalized deep water
 
solution for different ice thickness <math>0.1</math>, <math>1.0</math> and <math>10</math> metres. We find
 
that the positions of the roots (wavenumbers) of the non-dimensional
 
deep-water dispersion equation remain virtually unchanged as ice thickness is varied.
 
 
\begin{figure}[ptbh]
 
\begin{center}
 
\includegraphics[height=6.6206cm,width=7.8002cm]{scaleroots.eps}
 
\caption{Graph of the real and imaginary parts of the roots of the deep water dispersion equation agianst the non-dimensional radial frequency <math>\omega</math> plotted in a loglog scale. The thickness of the ice sheet is taken to be <math>0.1</math>, <math>1</math> and <math>10</math> metres.}
 
(fig:3-5)
 
\end{center}
 
\end{figure}
 
 
By closer observation of Fig.~((fig:3-5)), we find that up to
 
<math>\omega\approx1</math> all the roots are nearly identical for all the thickness of
 
the ice sheet considered. Although, the value of the roots other than
 
<math>q_{=T= }</math> vary slightly for higher frequencies, the position of the roots
 
stays qualitatively unchanged. Hence, we may conclude that our normalizing
 
scheme is valid for the range of ice thickness and frequency that are
 
geophysically relevant.
 
 
It is obvious that the same conclusion can be made for finite water depth,
 
since the depth dependent term, <math>\tanh\gamma H</math> in the dispersion equation is
 
dimensionless and independent of changes in ice-thickness.
 
Fig.~((fig:3-17)) shows the positions of the roots when the water depth
 
affects the response of the ice sheet, <math>H=0.2\pi</math> and <math>\pi</math>.
 
 
\begin{figure}[ptbh]
 
\begin{center}
 
\includegraphics[height=6.491cm,width=12.903cm]{shroots.eps}
 
\caption{Graphs of the real and imaginary parts of the roots of the shallow-water dispersion equation. (a) is when <math>H=0.2\pi</math> and (b) is when <math>H=\pi</math>.}
 
(fig:3-17)
 
\end{center}
 
\end{figure}
 
 
So far we have computed the deflection of the ice sheet with <math>m=0</math>. This is
 
justified because for typical ice sheets, it is reasonable to say that <math>m</math> is
 
smaller than <math>0.1</math> assuming that ice thickness is in the order of metres and
 
effective Young's modulus is in the order of <math>10^{9}</math> Pascals, and thus the
 
characteristic length <math>l_{=c= }</math> is typically larger than <math>10h^{3/4}</math>.
 
Furthermore <math>m</math> that is smaller than <math>0.1</math> has no noticeable effects on the
 
response of the ice sheet since the positions of the roots are virtually
 
unchanged as shown in Fig.((fig:3-18)). Omitting the mass density term,
 
the resulting dispersion equation for the deep water problem then becomes
 
truly independent of any physical parameters, which is the objective of the non-dimensionalization.
 
 
\begin{figure}[ptbh]
 
\begin{center}
 
\includegraphics[height=2.8037in,width=2.6481in]{scalrt.eps}
 
\caption{Positions of the roots of the non-dimensional dispersion equation when the non-dimensional mass density <math>m</math> is set zero and <math>0.1</math>. The radial frequency is ranging from <math>0.1</math> to <math>10</math>. }
 
(fig:3-18)
 
\end{center}
 
\end{figure}
 
 
====General scaling law of a floating ice sheet====
 
 
We have seen that the scaling scheme that was originally based on the static
 
loading of an infinite plate can be extended to dynamic problems. It is shown
 
by Fox [[Fox00]] that our scaling scheme can actually be extended to a
 
plate of finite size and the scaling law for the dynamics of an elastic
 
floating plate (ice sheet) can be found in a more general context. It can be
 
shown that the representation of the fundamental solution by a modal expansion
 
is independent of the shape of the plate, as long as its shape is reasonably
 
smooth, since the dispersion function, and hence the denominator of the
 
Fourier transform of <math>\bar{w}</math> is unaffected by the shape of the plate. This
 
subsection shows a brief description of the work done by C. Fox for the
 
article [[Fox00]], as it shows the relevance of the scaling regime in the
 
general physical situations.
 
 
Reviewing the method of solution of the infinite plate, we notice that the
 
geometry of the plate is unnecessary to obtain the dispersion function
 
<math>d\left(  \gamma,\omega\right)  </math>, since it is determined only by the plate equation.
 
 
We consider a smoothly shaped plate denoted by <math>\Omega</math>, as in the previous
 
chapter. Using the same notations in chapter 2, we are able to find a solution
 
of the Laplace's equation in <math>\mathcal{V}</math>. Then the Fourier transform of the
 
plate equation must integrated over the finite domain <math>\Omega</math>, which consists
 
of the terms to be determined by the values of <math>w</math> and its derivatives at the
 
boundary of the plate <math>\partial\Omega</math>.
 
<center><math>
 
\hat{\bar{w}}_{\Omega}\left(  \bar{\alpha},\bar{k}\right)  =\int_{\Omega}
 
\bar{w}\left(  \bar{x},\bar{y}\right)  e^{\mathrm{i}\left(  \bar
 
{\alpha}\bar{x}+\bar{k}\bar{y}\right)  }d\bar{x}d\bar{y}
 
</math></center>
 
Then, using the inverse Fourier transform to calculate the displacement,
 
<center><math>
 
\bar{w}\left(  \bar{x},\bar{y}\right)  =\frac{1}{4\pi^{2}}\int_{-\infty
 
}^{\infty}\int_{-\infty}^{\infty}\hat{\bar{w}}_{\Omega}\left(  \bar{\alpha
 
},\bar{k}\right)  e^{-\mathrm{i}\left(  \bar{\alpha}\bar{x}+\bar{k}
 
\bar{y}\right)  }d\bar{\alpha}d\bar{k}. (3-43)
 
</math></center>
 
Note that <math>w\left(  \bar{x},\bar{y}\right)  </math>\ calculated using Eqn.~((3-43)) is zero outside the ice sheet <math>\Omega</math>.
 
 
Using the Fourier transform defined above and Green's theorem, we find that the Fourier transform of the plate equation becomes
 
<center><math>
 
\left(  D\bar{\gamma}^{4}+\rho g-\bar{m}\bar{\omega}^{2}\right)  \hat{\bar{w}
 
}_{\Omega}\left(  \bar{\alpha},\bar{k}\right)  +\mathrm{i}\bar{\omega
 
}\hat{\bar{\phi}}_{\Omega}\left(  \bar{\alpha},\bar{k},\bar{z}\right)
 
=\bar{p}_{=a= }+b_{1}\left(  \bar{\alpha},\bar{k}\right)  (3-34)
 
</math></center>
 
where <math>b_{1}</math> is function of <math>\left(  \bar{\alpha},\bar{k}\right)  </math>. <math>b_{1}</math> is determined by the boundary values of <math>w</math> and its derivatives, which results Green's theorem. The Fourier transform of Laplace's equation in <math>\Omega</math> for <math>-\bar{H}<\bar{z}<0</math>\ becomes
 
<center><math>
 
\frac{\partial^{2}\hat{\bar{\phi}}_{\Omega}}{\partial z^{2}}-\bar{\gamma}
 
^{2}\hat{\bar{\phi}}_{\Omega}=b_{2}\left(  \bar{\alpha},\bar{k}\right)
 
(3-33)
 
</math></center>
 
where <math>b_{2}\left(  \bar{\alpha},\bar{k}\right)  </math> is again an inhomogeneous
 
term that arises in Green's theorem. We can solve Eqn.~((3-33)), which is
 
an ordinary differential equation, given the fixed ocean floor condition and
 
<math>b_{2}</math>. Thus, we find the relation between <math>\hat{\bar{\phi}}_{\Omega} </math> and
 
its <math>\bar{z}</math> derivative similar to that in the previous sections,
 
<center><math>
 
\frac{\partial\hat{\bar{\phi}}_{\Omega}}{\partial\bar{z}}\left(  \bar{\alpha
 
},\bar{k},0\right)  =\bar{\gamma}\tanh\left(  \bar{\gamma}\bar{H}\right)
 
\hat{\bar{\phi}}_{\Omega}\left(  \bar{\alpha},\bar{k},0\right)  +b_{3}\left(
 
\bar{\alpha},\bar{k}\right)
 
</math></center>
 
where <math>b_{3}</math> depends on the functions <math>b_{1}</math> and <math>b_{2}</math> in Eqn.~((3-34)) and Eqn.~((3-33)). Therefore the Fourier transform <math>\hat{\bar{w}}\left( \bar{\alpha},\bar{k}\right)  </math> has the same denominator, the dispersion function
 
<center><math>
 
\bar{d}\left(  \bar{\gamma},\bar{\omega}\right)  =D\bar{\gamma}^{4}+\rho
 
g-\bar{m}\bar{\omega}^{2}-\frac{\rho\bar{\omega}^{2}}{\bar{\gamma}\tanh
 
\bar{\gamma}\bar{H}}
 
</math></center>
 
as the infinite plate solution,
 
<center><math>
 
\hat{\bar{w}}\left(  \bar{\alpha},\bar{k}\right)  =\frac{b\left(  \bar{\alpha
 
},\bar{k}\right)  }{\bar{d}\left(  \bar{\gamma},\bar{\omega}\right)  }.
 
</math></center>
 
The inverse Fourier transform of <math>\hat{\bar{w}}\left(  \bar{\alpha},\bar
 
{k}\right)  </math>\ will then be dependent on the shape of the plate and the
 
boundary conditions given at the edge of the plate, which determine the
 
function <math>b\left(  \bar{\alpha},\bar{k}\right)  </math>. If we assume that the
 
numerator <math>b\left(  \bar{\alpha},\bar{k}\right)  </math> has no singularities, in
 
other words <math>\bar{w}</math> and its derivatives on the edge of the plate are
 
well-behaved functions, then the inverse Fourier transform of <math>\hat{\bar{w}
 
}\left(  \bar{\alpha},\bar{k}\right)  </math> must again be expressed by the mode
 
expansion over the roots of the dispersion equation, <math>K^{}</math>.
 
Because our scaling is based on the behavior of the roots of the dispersion
 
equation, the scaling can be applied to an ice sheet of general shape.
 
 
For simple shapes of the plate, for example a circular disk, semi-infinite
 
straight edged plate, or rectangular plate under line load parallel to the
 
edge of the plate, the function <math>b</math> will depend only on <math>\bar{\gamma}</math>, the
 
amplitude of <math>\left(  \bar{\alpha},\bar{k}\right)  </math>. Hence, the inverse
 
Fourier transform will again then depend only on the roots of the dispersion
 
equation, which can be effectively scaled using the characteristic length and
 
time. For an infinite plate, we have found that <math>b\equiv1</math>. We will see an
 
example of <math>b\neq1</math>, i.e., non-infinite plate, when the inverse Fourier
 
transform can be analytically calculated in the next chapter.
 
 
===Determining characteristic length from field measurements (sec:1)===
 
 
We consider how the characteristic length <math>l_{=c= }</math>, and hence the
 
effective Young's modulus, can be determined from flexural motion of the ice
 
sheet. As mentioned before, although the ice sheet is modeled with a constant
 
Young's modulus, the value of <math>E</math> varies through the ice sheet in reality
 
mainly due to the temperature gradient from top to bottom. As a result, we use
 
a constant ''effective''  Young's modulus as a substitute. A commonly used
 
value of the effective Young's modulus is <math>5\times10^{9}N
 
m  ^{-2}</math>. However, it is not obvious how one can obtain the actual Young's modulus and then the effective Young's modulus. We propose here a few possible methods of measuring the effective Young's modulus from field experiments using the theoretical results in the previous sections.
 
 
====Characteristic length====
 
 
Fig.~((fig:thumper)) shows a schematic drawing of how flexural waves in an ice sheet can be generated using a hydraulic "thumper" that can lift a block of ice up and down at a prescribed frequency. The amplitude of the oscillating force is calculated from the size of the ice block and the length of the stroke of the thumper.
 
 
\begin{figure}[ptbh]
 
\begin{center}
 
\includegraphics[height=3.2598cm,width=7.831cm]{thumper.eps}
 
\caption{Schematic of an experimented setup of a `thumper'. A block of ice is cut out and lifted up and down.}
 
(fig:thumper)
 
\end{center}
 
\end{figure}
 
 
We consider methods of calculating the characteristic length in two cases, when the water is shallow, i.e., <math>H=0.2\pi</math>, and when the water is deep, i.e., <math>H</math> larger than <math>2\pi</math>. We have seen that water depth larger than <math>2\pi</math> can be regarded as deep in the previous sections. Fig.~((fig:3-10)) and Fig.~((fig:3-11)) show the magnitude and phase of the normalized strain, which is the second derivative of the non-dimensional displacement function, <math>-w_{rr}\left(  r,\omega\right)  </math>, plotted as a function of non-dimensional radian frequency <math>\omega</math> and distance <math>r</math> from the point of forcing in the case when the water depth is <math>H=0.2\pi</math>. Then, the strain per Newton denoted by <math>\bar{S}\left(  \bar{r}\right)  </math> is calculated by
 
<center><math>\begin{matrix}
 
\bar{S}\left(  \bar{r}\right)  & =-\frac{h}{2}\bar{w}_{\bar{r}\bar{r}}\left(
 
\bar{r}\right)  =-\frac{h}{2\rho gl_{=c= }^{4}}w_{rr}\left(  r\right)
 
\\
 
& =-\frac{ih}{4\rho gl_{=c= }^{4}}\sum_{q\in K^{}}qR\left(
 
q\right)  \left(  q^{2}H_{2}^{\left(  1\right)  }\left(  qr\right)  -\frac
 
{q}{r}H_{1}^{\left(  1\right)  }\left(  qr\right)  \right)  . (3-66)
 
\end{matrix}</math></center>
 
We used the formula for the derivative of Hankel function (Abramowitz and
 
Stegun [[Stegun]]). The following graphs of the amplitude of the strain
 
shows the strain scaled by <math>h/\left(  2\rho gl_{=c= }^{=4= }\right)
 
<math>, hence </math>\left|  S\left(  r\right)  \right|  =\left|  w_{rr}\left(  r\right)
 
\right|  </math>.
 
 
\begin{figure}[ptb]
 
\begin{center}
 
\includegraphics[height=8.1012cm,width=8.7601cm]{strabs02.eps}
 
\caption{Magnitude of the normalized strain for shallow water <math>H=0.2\pi</math>, as  a function of non-dimensional radial forcing fequency (log-scale) and non-dimensional distance, <math>0.1\leq r\leq1.5</math>, from the forcing point (in reversed axis for better view).}
 
(fig:3-10)
 
\end{center}
 
\end{figure}
 
 
\begin{figure}[ptb]
 
\begin{center}
 
\includegraphics[height=8.1012cm,width=8.7601cm]{strph02.eps}
 
\caption{Phase of the strain for shallow water <math>H=0.2\pi</math>, as a function of non-dimensional radial forcing frequency (log-scale) and non-dimensional distance, <math>0.1\leq r\leq1.5</math>, from the forcing point.}
 
(fig:3-11)
 
\end{center}
 
\end{figure}
 
 
\begin{figure}[ptb]
 
\begin{center}
 
\includegraphics[height=8.1012cm,width=8.7601cm]{fig4.eps}
 
\caption{Magnitude of the normalized strain for deep water, as \ a function of non-dimensional radial forcing fequency (log-scale) and non-dimensional distance, <math>0.1\leq r\leq1.5</math>, from the forcing point (in reversed axis for better view).}
 
(fig:3-8)
 
\end{center}
 
\end{figure}
 
 
\begin{figure}[ptb]
 
\begin{center}
 
\includegraphics[height=8.1012cm,width=8.7601cm]{fig5.eps}
 
\caption{Phase of the strain for deep water, as a function of forcing frequency and distance, <math>0.1\leq r\leq1.5</math>, from the forcing point.}
 
(fig:3-9)
 
\end{center}
 
\end{figure}
 
 
Fig.~((fig:3-8)) and Fig.~((fig:3-9)) are equivalent plots for water
 
depth <math>H=2\pi</math>. Note that the plots of strain magnitude and phase have
 
reversed distance axes in order to show the structure of the curves.
 
Fig.~((fig:3-10)) and Fig.~((fig:3-8)) show that the magnitude of
 
strain at the near field, <math>r<1</math>, changes rapidly with <math>r</math>. Hence, because of
 
the length of the strain gauge itself and the physical size of the forcing
 
device, it is impossible to measure the strain at a point sufficiently
 
accurately to fit to the graph of magnitude. If instead the strain gauges are
 
placed near <math>r\gtrsim1</math>, where the magnitude of strain changes little with
 
respect to distance from the forcing then, by sweeping the frequency, we
 
should be able to find the frequency where the ''dip''  in magnitude occurs
 
and hence determine the characteristic frequency. A more robust measure that
 
does not require an accurate knowledge of the magnitude of forcing, is to use
 
the feature in Fig.~((fig:3-11)) and Fig.~((fig:3-9)) which show that
 
for <math>r\gtrsim 1</math> the phase of strain has a minimum value at <math>\omega\approx0.9</math> and the position of the minimum is an insensitive function of distance. Hence, robust measurements can be made at <math>r\gtrsim1</math> to find a frequency where the ''dip''  in phase occurs to determine characteristic time and length.
 
 
Fig.~((fig:3-7)) shows the magnitude and phase of the acceleration at the
 
point of forcing as a function of non-dimensional radial frequency <math>\omega</math>.
 
The acceleration is calculated from the formula of the displacement at the
 
point of forcing, Eqn.~((3-45)), <math>-\omega^{2}w\left(  0\right)  </math>. We
 
notice that the both magnitude and phase of the acceleration are monotonically
 
increasing functions of <math>\omega</math>. Hence, we may sweep frequency and measure
 
the response to match the theoretical results in Fig.~((fig:3-7)), to find
 
the characteristic length.
 
 
\begin{figure}[ptb]
 
\begin{center}
 
\includegraphics[height=8.1297cm,width=9.9112cm]{fig3.eps}
 
\caption{(a) Magnitude and (b) phase of normalized vertical acceleration at the point of forcing, as a function of non-dimensional frequency.}
 
(fig:3-7)
 
\end{center}
 
\end{figure}
 
 
It may be difficult to measure the magnitude of the acceleration, since it
 
requires an accurate knowledge of the weight of the ice block. Therefore, it
 
may be more robust to make use of the graphed phase of the acceleration. In
 
order to find the characteristic length, we sweep frequency to locate the
 
relative phase of <math>0.9</math> at which the shallow and deep water phases coincide
 
(at about <math>\omega=0.9</math>) so that water depth does not affect much. Let
 
<math>\bar{\omega}_{0}</math> denote such frequency, then we find <math>\bar{\omega}_{0}</math> at
 
which the measured relative phase is <math>0.9</math>, i.e.,
 
<center><math>
 
\arg\left(  \bar{w}\left(  0,\bar{\omega}_{0}\right)  \right)  =0.9.
 
</math></center>
 
From <math>\bar{\omega}_{0}</math> we are then able to find the characteristic length
 
using
 
<center><math>
 
l_{=c= }=\frac{0.9^{2}}{\bar{\omega}_{0}^{2}}g. (3-70)
 
</math></center>
 
 
Another possible measurement method of the characteristic length is to use the tilt of the ice sheet in far field calculated using an asymptotic expression
 
<center><math>
 
\frac{d}{dr}w\left(  r\right)  \sim-\frac{\mathrm{i}q_{=T= }
 
^{3/2}R_{=T= }}{\sqrt{2\pi r}}=\frac{-\mathrm{i}}{\sqrt{2\pi r}
 
}\frac{q_{=T= }^{7/2}}{5\omega^{2}-4q_{=T= }}=\frac{-\mathrm{i}
 
}{\sqrt{2\pi r}}\frac{q_{=T= }^{5/2}}{5q_{=T= }^{4}+1}.
 
</math></center>
 
We use the asymptotic formula of a Hankel function,
 
<center><math>
 
H_{\nu}^{\left(  1\right)  }\left(  r\right)  \sim\sqrt{\frac{2}{\pi r}}
 
\exp\left[  -\mathrm{i}\left(  r-\frac{\nu}{2}-\frac{\pi}{4}\right)
 
\right]  .
 
</math></center>
 
Again we try to locate the physical radial frequency <math>\bar{\omega}</math> at which
 
the maximum tilt occurs (at unit non-dimensional radial frequency <math>\omega=1</math>)
 
as seen in Fig.~((fig:3-27)). Then, using Eqn.~((3-70)), we find the
 
characteristic length.
 
 
\begin{figure}[ptb]
 
\begin{center}
 
\includegraphics[height=7.8683cm,width=9.8453cm]{tilt.eps}
 
\caption{Sacled amplitude of tilt as a function non-dimensional radial frequency <math>\omega=\bar{\omega}t_{=c= }</math> by lifted ice-weight (<math>10^{3}</math>kg) devided by <math>l_{=c= }^{3}</math> and square root of non-dimensional distance <math>\sqrt{r}</math>, <math>10^{3}=kg= /l_{=c= }^{3}=m= ^{3}/\sqrt{r}</math>.}
 
(fig:3-27)
 
\end{center}
 
\end{figure}
 
 
The methods of measuring the characteristic length proposed here assume that the forcing is localized only on the ice sheet. However, in practice it may be difficult to generate a pure surface-load on the ice because of the size of an ice-block which can be almost 2 metres tall. When the thumper is configured as in Fig.~((fig:thumper)), the effects of water (about 1 cubic metre) which is pumped in and out by the ice-block cannot be ignored. Therefore, if a thumper shown in Fig.~((fig:thumper)) is used to generate flexural motions a mathematical model that incorporates weight of an ice block and water-pumping action will be required to analyze experiment data sets. This is ongoing research.
 
 
===Summary===
 
 
We may divide the content of this chapter into two parts. One that has to do with mathematical calculations involving integration, complex valued functions, to special functions, and so on. Another that applies the analytical solutions to real geophysical studies of sea ice, e.g. for the effective scaling of the solution and in-situ measurement of Young's modulus.
 
 
The response of the ice sheet to a localized force is expressed by an  infinite series of natural wave modes of the ice sheet and the wavenumbers of the modes are the roots of the dispersion equation which is a relationship between the forcing frequency and the wavenumbers. The dispersion equation is extended to the complex plane by analytic continuation, which enables us to use the tools of contour integrals in the complex plane and commonly available tables of integral transformations. The process of deriving the formulae of the coefficients is tedious but a series of elementary calculations.
 
 
We are able to calculate the inverse transform analytically using the positions of the roots of the dispersion equation in the complex plane. When the water depth is finite, the solutions are expressed by infinite series of natural wave modes existing in an ice sheet. For infinite water depth case, the finite water depth solutions can be directly extended to derive a simple formula for the solution made up of five wave modes. Our derivation of the deep-water solution avoids the use of analyticity properties of the dispersion equation for deep-water problems. The formulae are simplified by non-dimensionalization of the equations and a consequent reduction of the number of physical parameters. We found that although the characteristic length <math>l_{=c= }</math> is originally based on the static solution, together with the characteristic time <math>t_{=c= }</math>, the dynamic system can also be effectively non-dimensionalized, so that the resulting non-dimensional solutions are insensitive to changes of physical parameters.
 
 
Knowing the theoretical conversion relationship between the scaled and actual solutions lets us find meaningful values of various parameters that can be interpreted to any physical scale using the characteristic length and time. For example, we found that non-dimensional water depth <math>H=2\pi</math> can be considered as deep regardless of the wavelength and frequency of surface waves, and that the maximum coupling between the forcing and the response of the ice sheet happens just below unit non-dimensional frequency <math>\omega=1</math>. Furthermore, our scaling scheme is shown to be applicable to more generally shaped ice sheet, hence the values of the non-dimensional parameters, such as <math>H=2\pi</math> and <math>\omega=1</math>, become important for non-infinite ice sheets.
 
 
Section 3.6 proposes several possible methods for how the theoretical solution can be used to obtain an effective Young's modulus from an actual field experimental data set. Experimental data show that existing strain gauges can be used to determine characteristic length using the measurement schemes introduced in section 3.6.
 
 
All plots of displacement are computed using the software package MatLab on an Intel Pentium III PC and all the graphs shown in this chapter are the results of the direct implementation of the formulae given here. The number of roots used for the computation is proportional to the water depth,
 
<center><math>
 
=number of roots= \approx\left(  \frac{H}{2\pi}+2\right)  \times10,
 
</math></center>
 
which is about <math>20</math> to <math>100</math> for each given frequency. The number of roots given by the formula above may be larger than the minimum number required, but finding the roots is computationally inexpensive. The special functions other than Struve function and <math>Ci </math> and <math>si </math> are computed by the built-in functions of MatLab. Struve function and <math>Ci </math> and <math>si </math> are computed using the power expansion given in appendix A.
 
 
The idea from the dispersion equation to the series expansion Eqn.~((3-50)) was given to me by my supervisor C. Fox. The computation of the Green's function (except the roots of the dispersion equation) for the finite and the infinite water depth cases are done by me.
 
 
  
  
 
[[Category:Floating Elastic Plate]]
 
[[Category:Floating Elastic Plate]]

Latest revision as of 16:40, 8 December 2009


This is a special version of the free-surface Green function which applied when the Floating Elastic Plate boundary condition applies at the free-surface. It reduced to the Free-Surface Green Function in the limit as the plate terms tend to zero and to the Green function for an infinite plate in the limit as the water terms tend to zero.

Two Dimensions

Green function with singularity on the surface

For the case of a floating thin plate we almost always make the assumption of shallow draft and the Green function satisfies the following equations The Green function [math]\displaystyle{ G(x,x^{\prime},z) }[/math] representing outgoing waves as [math]\displaystyle{ |x|\rightarrow \infty }[/math] satisfies

[math]\displaystyle{ \nabla^2 G = 0, -H\lt z\lt 0, }[/math]
[math]\displaystyle{ \frac{\partial G}{\partial z} =0, z=-H, }[/math]
[math]\displaystyle{ {\left( \beta \frac{\partial^4}{\partial x^4} - \gamma\alpha + 1\right)\frac{\partial G}{\partial z} - \alpha G} = \delta(x-x^{\prime}), z=0, }[/math]

(note we are only considering singularities at the free surface). It is important to realise that these equations are based on the Frequency Domain Problem and that the exact form of the Green function is dependent on whether we have [math]\displaystyle{ e^{i\omega t} }[/math] or [math]\displaystyle{ e^{-i\omega t} }[/math] dependence. In what follows we assume [math]\displaystyle{ e^{i\omega t} }[/math] dependence.

This problem can be solved (a solution is given below and also in Evans and Porter 2006) to give

[math]\displaystyle{ G(x,x^{\prime},z) = -\sum_{n=-2}^\infty\frac{\sin{(k_n H)}\cos{(k_n(z+H))}}{2\alpha C(k_n)}e^{-k_n|x-x^{\prime}|}, }[/math]

where

[math]\displaystyle{ C(k_n)=\frac{1}{2}\left(h - \frac{(5\beta k_n^4 + 1 - \alpha\gamma)\sin^2{(k_n H)}}{\alpha}\right), }[/math]

and [math]\displaystyle{ k_n }[/math] are the solutions of the Dispersion Relation for a Floating Elastic Plate,

[math]\displaystyle{ \beta k^5 \sin(kH) - k \left(1 - \alpha \gamma \right) \sin(kH) = -\alpha \cos(kH) \, }[/math]

with [math]\displaystyle{ n=-1,-2 }[/math] corresponding to the complex solutions with positive real part, [math]\displaystyle{ n=0 }[/math] corresponding to the imaginary solution with positive imaginary part, since we have assumed [math]\displaystyle{ e^{i\omega t} }[/math] dependence (it would be of negative imaginary part if we had chosen [math]\displaystyle{ e^{-i\omega t} }[/math] dependence) [math]\displaystyle{ n\gt 0 }[/math] corresponding to the real solutions with positive real part. We can also write the Green function as

[math]\displaystyle{ G\left( r,z\right) = -\omega\sum_{n=-2}^{\infty} \frac{R\left( q_n\right) }{q_n\sinh q_nH}\exp\left( \mathrm{i}q_nr\right) \cosh q_n\left( z+H\right) . }[/math]

where [math]\displaystyle{ r=|x-x^{\prime}| }[/math], [math]\displaystyle{ q_n=ik_n }[/math] and we have used a different non-dimensionalisation so that the Dispersion Relation for a Floating Elastic Plate is written as

[math]\displaystyle{ -q^5 \sinh(kH) - k \left(1 - m\omega^2 \right) \sinh(kH) = -\omega^2 \cosh(kH) \, }[/math]

and

[math]\displaystyle{ R\left( q_n\right) =\frac{\omega^{2}q_n}{\omega^{2}\left( 5q_n^{4}+u\right) +H\left[ \left( q_n^{5}+uq\right) ^{2}-\omega^{4}\right] }. }[/math]

Green function with singularity in the water

We can define the Green function [math]\displaystyle{ G(x,x^{\prime},z) }[/math], representing outgoing waves as [math]\displaystyle{ |x|\rightarrow \infty }[/math], which satisfies

[math]\displaystyle{ \nabla^2 G = \delta(x-x^{\prime})\delta(z-z^{\prime}), -H\lt z\lt 0, }[/math]
[math]\displaystyle{ \frac{\partial G}{\partial z} =0, z=-H, }[/math]
[math]\displaystyle{ {\left( \beta \frac{\partial^4}{\partial x^4} - \gamma\alpha + 1\right)\frac{\partial G}{\partial z} - \alpha G} = 0, z=0. }[/math]

Note that this Green function reduces the Free-Surface Green Function in the limit as [math]\displaystyle{ \beta }[/math] and [math]\displaystyle{ \gamma }[/math] tend to zero.

This problem can be solved to give Evans and Porter 2003

[math]\displaystyle{ G(x,x^{\prime},z) = -\sum_{n=-2}^\infty\frac{\cos{(k_n(z^{\prime}+H))}\cos{(k_n(z+H))}}{2k_n C(k_n)}e^{-k_n|x-x^{\prime}|}, }[/math]

where [math]\displaystyle{ C(k_n) }[/math] and [math]\displaystyle{ k_n }[/math] are as before.

Relationship between the two Green functions

If we denote the Green function with the singularity at the free-surface by [math]\displaystyle{ G_s }[/math] and the Green function with the singularity in the water by [math]\displaystyle{ G_w }[/math]then we have the following relationship between them

[math]\displaystyle{ G_s = -\frac{1}{\alpha}\left. \frac{\partial G_w}{\partial z^{\prime}} \right|_{z^{\prime} =0} }[/math]

Three Dimensions

[math]\displaystyle{ G\left( r,z\right) =-\frac{\omega}{2}\sum_{n=-2}^{\infty} \frac{R\left( q_n\right) }{\sinh q_nH}H_{0}^{\left( 1\right) }\left( q_n r\right) \cosh q_n\left( z+H\right) }[/math]

Derivation of the Green function

We present here a derivation of the Green function for both the two and three dimensional problem. The devivation uses the Mittag-Leffler Expansion for the Floating Elastic Plate Dispersion Equation. This derivation is based on the PhD thesis of Hyuck Chung.

Previous Work

In the following sections, the deflection of a plate is expressed by the infinite number of modes that exist in the vibrating plate floating on the water, instead of a finite number of modes that represent the free oscillation of an elastic plate without restraints. The fact that waves in an incompressible fluid can be expressed by an infinite series of natural modes is shown by John 1949 and John 1950 and discussed in Cylindrical Eigenfunction Expansion Each mode being a Hankel function of the first kind in the horizontal (in the three-dimensional case) and a hyperbolic function in the vertical. Kheisin 1967 chapter IV, studied the same problem solved in this chapter and derived the same dispersion equation for the physical variables. He studied the properties of the inverse transform only for the simple shallow water and static load cases, which are approximated versions of the full solutions which is presented here.

Mathematical model

In terms of non-dimensional variables the Floating Elastic Plate equation becomes

[math]\displaystyle{ \left[ \nabla^{4}-m\omega^{2}+1\right] w+\mathrm{i}\omega \phi=p_n. }[/math]

where [math]\displaystyle{ p_n }[/math] is the pressure and [math]\displaystyle{ n }[/math] is the dimension. We assume that [math]\displaystyle{ p_2 = \delta(x) }[/math] and that [math]\displaystyle{ p_3 = \delta(x) \delta(y) }[/math]. [math]\displaystyle{ w }[/math] is the displacement of the surface which is related to the velocity potential by the equation [math]\displaystyle{ \frac{\partial\phi}{\partial z} = \frac{\partial w}{\partial t} }[/math] We also have Laplace's equation for the velocity potential and the bottom condition

[math]\displaystyle{ \nabla^{2}\phi=0 }[/math]

for [math]\displaystyle{ -\infty\lt x,y\lt \infty, }[/math] [math]\displaystyle{ -H\lt z\lt 0 }[/math] and

[math]\displaystyle{ \phi_{z}=0\quad \mathrm{at} \quad \;z=-H. }[/math]

The system of equations from together with the Sommerfeld Radiation Condition form the BVP, which we will solve for [math]\displaystyle{ w }[/math] and [math]\displaystyle{ \phi }[/math] for a given [math]\displaystyle{ \omega }[/math].

Spatial Fourier transform

We solve the system of equations using the Fourier Transform in [math]\displaystyle{ \left( x,y\right) }[/math]-plane for three-dimensions and with a Fourier Transform in [math]\displaystyle{ x }[/math] for two dimensions. We choose the Fourier transform with respect to [math]\displaystyle{ x }[/math] and [math]\displaystyle{ y }[/math] defined as

[math]\displaystyle{ \hat{\phi}\left( \alpha,k,z\right) =\int_{-\infty}^{\infty}\int_{-\infty }^{\infty}\phi\left( x,y,z\right) e^{\mathrm{i}\left( \alpha x+ky\right) }\mathrm{d}x\mathrm{d}y }[/math]

and the inverse Fourier Transform defined as

[math]\displaystyle{ \phi\left( x,y,z\right) =\frac{1}{4\pi^{2}}\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\hat{\phi}\left( \alpha,k,z\right) e^{-\mathrm{i}\left( \alpha x+ky\right) }\mathrm{d}\alpha \mathrm{d}k. }[/math]

For the one-dimensional case, the definitions are

[math]\displaystyle{ \hat{\phi}\left( \alpha,z\right) =\int_{-\infty}^{\infty}\phi\left( x,z\right) e^{\mathrm{i}\alpha x}\mathrm{d}x, }[/math]
[math]\displaystyle{ \phi\left( x,z\right) =\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat{\phi }\left( \alpha,z\right) e^{-\mathrm{i}\alpha x}\mathrm{d}\alpha. }[/math]

We denote the spatial Fourier Transform by using a hat over [math]\displaystyle{ w }[/math] and [math]\displaystyle{ \phi }[/math].

The Fourier Transform of both sides of the Laplace's equation becomes an ordinary differential equation with respect to [math]\displaystyle{ z }[/math],

[math]\displaystyle{ \frac{\partial^{2}\hat{\phi}}{\partial z^{2}}\left( \alpha,k,z\right) -\left( \alpha^{2}+k^{2}\right) \hat{\phi}\left( \alpha,k,z\right) =0 }[/math]

with a solution

[math]\displaystyle{ \hat{\phi}\left( \alpha,k,z\right) =A\left( \gamma\right) e^{\gamma z}+B\left( \gamma\right) e^{-\gamma z} }[/math]

where [math]\displaystyle{ \gamma^{2}=\alpha^{2}+k^{2} }[/math] and [math]\displaystyle{ A }[/math], [math]\displaystyle{ B }[/math] are unknown coefficients to be determined by the boundary condition. We find that [math]\displaystyle{ \hat{\phi} }[/math] is a function only of the magnitude of the Fourier transform variables, [math]\displaystyle{ \gamma=\left\| \mathbf{\gamma}\right\| =\left| \left| \left( \alpha,k\right) \right| \right| }[/math], hence we may now denote [math]\displaystyle{ \hat{\phi }\left( \alpha,k,z\right) }[/math] by [math]\displaystyle{ \hat{\phi}\left( \gamma,z\right) }[/math]. We can reach the same conclusion using the fact that [math]\displaystyle{ w\left( x,y\right) }[/math] and [math]\displaystyle{ \phi\left( x,y,z\right) }[/math] are functions of [math]\displaystyle{ r=\sqrt{x^{2}+y^{2}} }[/math] thus the Fourier transforms must be functions of [math]\displaystyle{ \gamma=\sqrt{\alpha^{2}+k^{2}} }[/math].

We find the unknown coefficients [math]\displaystyle{ A }[/math] and [math]\displaystyle{ B }[/math] from the Fourier transformed ocean floor condition that [math]\displaystyle{ \left. \hat{\phi}_{z}\right| _{z=-H}=0 }[/math] to be [math]\displaystyle{ A\left( \gamma\right) =Ce^{\gamma H} }[/math], and [math]\displaystyle{ B\left( \gamma\right) =Ce^{-\gamma H} }[/math]. Thus, we obtain the depth dependence of the Fourier transform of the potential

[math]\displaystyle{ \hat{\phi}\left( \gamma,z\right) =\hat{\phi}\left( \gamma,0\right) \frac{\cosh\gamma\left( z+H\right) }{\cosh\gamma H}. }[/math]

At the surface, [math]\displaystyle{ z=0 }[/math], differentiating both sides with respect to [math]\displaystyle{ z }[/math] the vertical component of the velocity is

[math]\displaystyle{ \hat{\phi}_{z}\left( \gamma,0\right) =\hat{\phi}\left( \gamma,0\right) \gamma\tanh\gamma H. }[/math]

Using this relationship to substitute for [math]\displaystyle{ \hat{\phi}_{z} }[/math] in the non-dimensional form of the kinematic condition, we find that

[math]\displaystyle{ \hat{\phi}\left( \gamma,0\right) =\frac{\mathrm{i}\omega\hat {w}\left( \gamma\right) }{\gamma\tanh\gamma H}. }[/math]

Thus, once we find [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] then we can find [math]\displaystyle{ \hat{\phi}\left( \gamma,z\right) }[/math].

The Fourier transform of the plate equation is also an algebraic equation in the parameter [math]\displaystyle{ \gamma }[/math]

[math]\displaystyle{ \left( \gamma^{4}+1-m\omega^{2}\right) \hat{w}+\mathrm{i}\omega \hat{\phi}=1. }[/math]

Hence, we have the single algebraic equation for each spatial Fourier component of [math]\displaystyle{ w }[/math]

[math]\displaystyle{ \left( \gamma^{4}+1-m\omega^{2}-\frac{\omega^{2}}{\gamma\tanh\gamma H}\right) \hat{w}=1. }[/math]

Notice that we have used the fact that the Fourier transform of the delta function is [math]\displaystyle{ 1 }[/math].

We find that the spatial Fourier transform of the displacement of the ice sheet for the localized forcing, both point and line, is

[math]\displaystyle{ \hat{w}\left( \gamma\right) =\frac{1}{d\left( \gamma,\omega\right) } }[/math]

where

[math]\displaystyle{ d\left( \gamma,\omega\right) =\gamma^{4}+1-m\omega^{2}-\frac{\omega^{2} }{\gamma\tanh\gamma H}. }[/math]

Where [math]\displaystyle{ d\left( \gamma,\omega\right) }[/math] in the Dispersion Relation for a Floating Elastic Plate function, and the associated algebraic equation [math]\displaystyle{ d\left( \gamma ,\omega\right) =0 }[/math] the dispersion equation for waves propagating through an ice sheet. This dispersion relation (and the Fourier transform of it) was derived by Kheisin 1967 chapter IV, and Fox and Squire 1994.

Inverse Fourier Transform

Our task now is to derive the inverse Fourier transform). We notice that the roots of the Dispersion Relation for a Floating Elastic Plate for a fixed [math]\displaystyle{ \omega }[/math], are the poles of the function [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math], which are necessary for calculation of integrals involved in the inverse Fourier transform. There are a pair of two real roots [math]\displaystyle{ \pm q_0 }[/math], an infinite number of pure imaginary roots [math]\displaystyle{ \pm q_{n},\,n\in\mathbb{N} }[/math], ([math]\displaystyle{ \mathrm{Im} q_{n}\gt 0 }[/math]) plus four complex roots (in all physical situations) occur at plus and minus complex-conjugate pairs [math]\displaystyle{ \pm q_{-1} }[/math] and [math]\displaystyle{ \pm q_{-2} }[/math] with [math]\displaystyle{ \mathrm{Im} \left( q_{-1,-2}\right) \gt 0 }[/math]).

We note that the dispersion equation represents a relationship between the spatial wavenumbers [math]\displaystyle{ \gamma }[/math] and the radial frequency [math]\displaystyle{ \omega }[/math], which is how the name "dispersion" came about.

We may also solve for the velocity potential at the surface of the water

[math]\displaystyle{ \hat{\phi}\left( \gamma,0\right) =\frac{\mathrm{i}\omega w} {\gamma\tanh\gamma H}=\frac{\mathrm{i}\omega}{\gamma\tanh\left( \gamma H\right) d\left( \gamma,\omega\right) }, (3-10) }[/math]

which is also a function of [math]\displaystyle{ \gamma }[/math] only and has the same poles as [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] does.

We calculate the displacement [math]\displaystyle{ w\left( x,y\right) }[/math] using the inverse Fourier transform of [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math]. Since [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] is radially symmetric, the inverse transform may be written (Bracewell 1965)

[math]\displaystyle{ w\left( r\right) =\frac{1}{2\pi}\int_{0}^{\infty}\hat{w}\left( \gamma\right) \gamma J_{0}\left( \gamma r\right) \mathrm{d}\gamma }[/math]

for the three dimensional problem, where [math]\displaystyle{ J_{0} }[/math] is Bessel function and [math]\displaystyle{ r }[/math] is the distance from the point of forcing. For the two-dimensional problem the inverse Fourier transform of [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] in the [math]\displaystyle{ x }[/math]-axis and since [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] is an even function, this is

[math]\displaystyle{ w\left( r\right) =\frac{1}{\pi}\int_{0}^{\infty}\hat{w}\left( \gamma\right) \cos\left( \gamma r\right) \mathrm{d}\gamma }[/math]

where again [math]\displaystyle{ r=\left| x\right| }[/math]. Note that the factors [math]\displaystyle{ 1/2\pi }[/math] and [math]\displaystyle{ 1/\pi }[/math] result from the form of the Fourier transform in two and one dimensional spaces.

The integrals are calculated using the Mittag-Leffler Expansion for the Floating Elastic Plate Dispersion Relation where we show that

[math]\displaystyle{ \hat{w}\left( \gamma\right) =\sum_{n=-2}^{\infty}\frac{2q_nR\left( q_n\right) }{\gamma^{2}-q_n^{2}} }[/math]

where [math]\displaystyle{ R\left( q_n\right) }[/math] is the residue of [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] at [math]\displaystyle{ \gamma=q_n }[/math] given by

[math]\displaystyle{ R\left( q_n\right) =\frac{\omega^{2}q_n}{\omega^{2}\left( 5q_n^{4}+u\right) +H\left[ \left( q_n^{5}+uq_n\right) ^{2}-\omega^{4}\right] }. }[/math]

.

Using the identities Abramowitz and Stegun 1964 formula 11.4.44 with [math]\displaystyle{ \nu=0 }[/math], [math]\displaystyle{ \mu=0 }[/math], [math]\displaystyle{ z=-\mathrm{i}q }[/math] and [math]\displaystyle{ a=r }[/math])

[math]\displaystyle{ \int_{0}^{\infty}\frac{\gamma}{\gamma^{2}-q^{2}}J_{0}\left( \gamma r\right) \mathrm{d}\gamma=K_{0}\left( -\mathrm{i}qr\right) }[/math]

for [math]\displaystyle{ Im q\gt 0 }[/math], [math]\displaystyle{ r\gt 0 }[/math], where [math]\displaystyle{ K_{0} }[/math] is a modified Bessel function and an identity between the modified Bessel function.

We can calculate the integral of the inverse Fourier transform of [math]\displaystyle{ \hat{w} }[/math] and, using the identity between [math]\displaystyle{ K_{0} }[/math] and Hankel function of the first kind (Abramowitz and Stegun 1964 formula 9.6.4),

[math]\displaystyle{ K_{0}\left( \zeta\right) =\frac{\mathrm{i}\pi}{2}H_{0}^{\left( 1\right) }\left( \mathrm{i}\zeta\right) }[/math]

for [math]\displaystyle{ Re \zeta\geq0 }[/math], the displacement for three-dimensions may be written in the equivalent forms

[math]\displaystyle{ w\left( r\right) =\frac{1}{\pi}\sum_{n=-2}^{\infty} q_n R\left( q_n\right) K_{0}\left( -\mathrm{i}q_n r\right) }[/math]
[math]\displaystyle{ =\frac{\mathrm{i}}{2}\sum_{n=-2}^{\infty}q_n R\left( q_n\right) H_{0}^{\left( 1\right) }\left( q_n r\right). }[/math]

In two-dimensions the identity (Erdelyi et. al.1954 formula 1.2 (11) with [math]\displaystyle{ x=\gamma }[/math], [math]\displaystyle{ y=r }[/math], and [math]\displaystyle{ a=-\mathrm{i}q }[/math])

[math]\displaystyle{ \int_{0}^{\infty}\frac{\cos\left( \gamma r\right) }{\gamma^{2}-q^{2}} \mathrm{d}\gamma=-\frac{\pi}{\mathrm{i}2q}\exp\left( \mathrm{i}qr\right) }[/math]

for [math]\displaystyle{ Im q\gt 0 }[/math], [math]\displaystyle{ r\gt 0 }[/math] can be used. This gives us

[math]\displaystyle{ w\left( r\right) =\mathrm{i}\sum_{n=-2}^{\infty} R\left( q\right) \exp\left( \mathrm{i}qr\right). }[/math]

Solution for the Velocity Potential

The velocity potential in the water can be found in two-dimensions

[math]\displaystyle{ \phi\left( r,z\right) =-\frac{\omega}{2}\sum_{n=-2}^{\infty} \frac{R\left( q_n\right) }{\sinh q_nH}H_{0}^{\left( 1\right) }\left( q_nr\right) \cosh q_n\left( z+H\right) }[/math]

and for three-dimensions

[math]\displaystyle{ \phi\left( r,z\right) =-\omega\sum_{n=-2}^{\infty} \frac{R\left( q_n\right) }{q_n\sinh q_nH}\exp\left( \mathrm{i}q_nr\right) \cosh q_n\left( z+H\right). }[/math]