Free-Surface Green Function for a 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
Two Dimensions
The Green function [math]\displaystyle{ G(x,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), z=0, }[/math]
(note we are only considering singularities at the free surface).
This problem can be solved to give
where
and [math]\displaystyle{ k_n }[/math] are the solutions of the Dispersion Relation for a Floating Elastic Plate, 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 and [math]\displaystyle{ n\gt 0 }[/math] corresponding to the real solutions with positive real part. We can also write the Green function as
where [math]\displaystyle{ q_n=ik_n }[/math] and
As each pole [math]\displaystyle{ q_n }[/math] is a root of the dispersion equation, we may substitute [math]\displaystyle{ \tanh q_nH=\omega^{2}/\left( q_n^{5}+uq_n\right) }[/math], where for brevity we have defined [math]\displaystyle{ u=\left( 1-m\omega^{2}\right) }[/math]. The residue may then be given as the rational function of the pole
Three Dimensions
Derivation of the Green function
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
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
whose solutions satisfy the radial form of the (non-dimensional) thin plate equation
The hydrodynamic effects that appear as [math]\displaystyle{ \phi_{t} }[/math] in the thin plate equation force us to solve Laplace's equation, and thus the application of the Fourier transform is not as straightforward as in the two examples cited (Sneddon Sneddon45 and Wyman Wyman50). Response of an infinite thin elastic plate under fluid loading is studied in a series of papers by Crighton Crighton72, Crighton77, Crighton79. 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 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 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 characteristic length, and thus the effective Young's modulus of a fast ice sheet using the non-dimensional solutions. In Fox, et al. Fox01, characteristic length and effective Young's modulus are estimated from actual experimental data sets collected by Fox at\ McMurdo Sound, Antarctica.
Fundamental solution
In this section, we derive a fundamental solution for a localized forcing with a single frequency applied on the ice sheet that is occupying the whole surface of the ocean. The fundamental solution for finite-depth water 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]\displaystyle{ w\left( x,y\right) }[/math], and then perform analytical integrations of the inverse Fourier transform.
Mathematical model
In terms of non-dimensional variables the Floating Elastic Plate equation becomes
Laplace's equation for the velocity potential and the bottom condition ((2-12)) remain the same, i.e.,
for [math]\displaystyle{ -\infty\lt x,y\lt \infty, }[/math] [math]\displaystyle{ -H\lt z\lt 0 }[/math] and
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
and the inverse Fourier transform defined as
For the one-dimensional case, the definitions are
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 ((3-2)) becomes an ordinary differential equation with respect to [math]\displaystyle{ z }[/math],
with a solution
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
At the surface, [math]\displaystyle{ z=0 }[/math], differentiating both sides of Eqn.~((3-5)) with respect to [math]\displaystyle{ z }[/math],\ the vertical component of the velocity is
Using this relationship to substitute for [math]\displaystyle{ \hat{\phi}_{z} }[/math] in the non-dimensional form of the kinematic condition, we find that
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] using Eqn.~((3-27)) and Eqn.~((3-6)).
The Fourier transform of the plate equation ((3-1)) is also an algebraic equation in the parameter [math]\displaystyle{ \gamma }[/math]
Hence, substituting Eqn.~((3-6)), we have the single algebraic equation for each spatial Fourier component of [math]\displaystyle{ w }[/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
where
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 was derived by Kheisin (Kheisin67 chapter IV), Fox and Squire fox3.
\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]\displaystyle{ \pm q_{=T= } }[/math] and on the left the imaginary roots [math]\displaystyle{ iq_{1},iq_{2},iq_{3},.... }[/math]. Dotted curves are [math]\displaystyle{ \tanh }[/math] function on the left and [math]\displaystyle{ \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]\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. 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]\displaystyle{ d\left( \gamma,\omega\right) }[/math] is an even function of [math]\displaystyle{ \gamma }[/math], and hence if [math]\displaystyle{ q }[/math] is a root then so is [math]\displaystyle{ -q }[/math]. Fig.~((fig:3-12)) shows how the real and imaginary roots can be found from the curves of the [math]\displaystyle{ \tan }[/math] and [math]\displaystyle{ \tanh }[/math] functions intersecting the curves
We notice that from observing the curves in Fig.~((fig:3-12)), there are a pair of two real roots [math]\displaystyle{ \pm q_{=T= } }[/math] ([math]\displaystyle{ q_{=T= }\gt 0 }[/math]) and an infinite number of pure imaginary roots [math]\displaystyle{ \left\{ \pm\mathrm{i}q_{n}\right\} _{n=1,2,...} }[/math], ([math]\displaystyle{ q_{n}\gt 0 }[/math]). We notice that for typical cases [math]\displaystyle{ m\omega^{2} \leq1 }[/math], an imaginary root [math]\displaystyle{ \mathrm{i}q_{n} }[/math] is found in each interval [math]\displaystyle{ \left( n-\frac{1}{2}\right) \pi\lt \gamma_{n}H\lt n\pi }[/math] as shown in Fig.~((fig:3-12)). For fixed [math]\displaystyle{ \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]\displaystyle{ \pm q_{=D= } }[/math] and [math]\displaystyle{ \pm q_{=D= }^{} }[/math] ([math]\displaystyle{ Re \left( q_{=D= }\right) \gt 0 }[/math] and [math]\displaystyle{ Im \left( q_{=D= }\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
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.
The Fourier transform of the differential equations can be interpreted as wave-like forcing of the ice sheet, i.e.,
where [math]\displaystyle{ \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]\displaystyle{ p_{=a= } }[/math] can be written as a function of [math]\displaystyle{ r }[/math],
where [math]\displaystyle{ \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]\displaystyle{ w\left( x,y\right) }[/math] by performing the two dimensional inverse Fourier transform of [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] in Eqn.~((3-8)) and, since [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] is radially symmetric, the inverse transform may be written (Bracewell Bracewell65)
where [math]\displaystyle{ J_{0} }[/math] is Bessel function and [math]\displaystyle{ r }[/math] is the distance from the point of forcing. The response to line forcing is given by 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
where again [math]\displaystyle{ r=\left| x\right| }[/math] is the distance from the line of forcing. 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 that we defined in Eqn.~((3-22)) and Eqn.~((3-36)).
The integrals of Eqn.~((3-11)) and Eqn.~((3-12)) are calculated using the singularities of [math]\displaystyle{ \hat{w} }[/math], i.e., the roots of the dispersion equation. First, since [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] is an even fractional function that equals zero when [math]\displaystyle{ \gamma=0 }[/math] and is bounded in the whole plane except in regions around its poles, we find that [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] can be expressed using the Mittag-Leffler expansion in section 3.2. Pairing each pole [math]\displaystyle{ q }[/math] with its negative counterpart [math]\displaystyle{ -q }[/math], gives
where [math]\displaystyle{ R\left( q\right) }[/math] is the residue of [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] at [math]\displaystyle{ \gamma=q }[/math]. We denoted the set of roots of the dispersion equation with positive imaginary part and a positive real root by [math]\displaystyle{ K^{} }[/math]. Thus, this set is [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ q\in K^{} }[/math].
\begin{figure}[tbh] \begin{center} \includegraphics[height=5.0918cm,width=7.5037cm]{polyres.eps} \caption{Graphs of [math]\displaystyle{ \left| R\left( \mathrm{i}q\right) \right| }[/math] as function of real number [math]\displaystyle{ q }[/math], with [math]\displaystyle{ \tanh }[/math]-functions (solid line) and polynomial expression (dashed line).}
(fig:3-31)
\end{center} \end{figure}
The residues [math]\displaystyle{ R\left( q\right) }[/math] can be calculated using the usual formula. Since each of the poles of [math]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] is simple, the residue [math]\displaystyle{ R\left( q\right) }[/math] at a pole [math]\displaystyle{ q }[/math] can be found using the expression
As each pole [math]\displaystyle{ q }[/math] is a root of the dispersion equation, we may substitute [math]\displaystyle{ \tanh qH=\omega^{2}/\left( q^{5}+uq\right) }[/math], where for brevity we have defined [math]\displaystyle{ u=\left( 1-m\omega^{2}\right) }[/math]. The residue may then be given as the rational function of the pole
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]\displaystyle{ q_{n}H }[/math] tends to [math]\displaystyle{ n\pi }[/math] as [math]\displaystyle{ n }[/math] becomes large, which makes [math]\displaystyle{ \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]\displaystyle{ \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]\displaystyle{ \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]\displaystyle{ \nu=0 }[/math], [math]\displaystyle{ \mu=0 }[/math], [math]\displaystyle{ z=-\mathrm{i}q }[/math] and [math]\displaystyle{ a=r }[/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 notice that [math]\displaystyle{ q_{=T= } }[/math] term of Eqn.~((3-50)) may pose a problem. However, considering that [math]\displaystyle{ 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]\displaystyle{ \mathrm{i}\beta\omega }[/math], [math]\displaystyle{ \beta\gt 0 }[/math], representing damping, to the dispersion equation
Every element in [math]\displaystyle{ K^{} }[/math] has positive imaginary part. Addition of the damping term proportional to [math]\displaystyle{ w_{t} }[/math] ensures that the inverse transform satisfies the radiation condition. We omit the damping term here for the sake of algebraic simplicity.
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 Stegun 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 point forcing may be written in the equivalent forms
where the subscript P of [math]\displaystyle{ w_{=P= } }[/math] indicates the response to a point load.
The identity (Erdelyi bateman formula 1.2 (11) with [math]\displaystyle{ x=\gamma }[/math], [math]\displaystyle{ y=r }[/math], and [math]\displaystyle{ a=-\mathrm{i}q }[/math])
for [math]\displaystyle{ Im q\gt 0 }[/math], [math]\displaystyle{ r\gt 0 }[/math]. The above integration gives the surface displacement for line forcing
where the subscript L of [math]\displaystyle{ 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]\displaystyle{ H }[/math] is found by first finding the roots in the upper-half plane, [math]\displaystyle{ 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
We used the identities [math]\displaystyle{ -\mathrm{i}\left( -q_{=D= }\right) =\left( -\mathrm{i}q_{=D= }\right) ^{}\lt math\gt , }[/math]R\left( q^{ }\right) =\left( R\left( q\right) \right) ^{}[math]\displaystyle{ 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]\displaystyle{ K^{} }[/math] by [math]\displaystyle{ R_{=T= }=R\left( q_{=T= }\right) }[/math], [math]\displaystyle{ R_{=D= }=R\left( q_{=D= }\right) }[/math], and [math]\displaystyle{ R_{n}=R\left( \mathrm{i}q_{n}\right) \lt math\gt , }[/math]n\in\mathbf{N}</math>, respectively. Note that [math]\displaystyle{ R\left( -\gamma_{=D= }^{}\right) \rightarrow-R_{=D= }^{} }[/math] as [math]\displaystyle{ \beta\rightarrow0 }[/math]. Since, from Eqn.~((3-17)), [math]\displaystyle{ \left| qR\left( q\right) \right| \propto\left| q\right| ^{-8}\lt math\gt 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]\displaystyle{ H }[/math], the fewer number of roots we need to achieve the desired accuracy, since for a given [math]\displaystyle{ n }[/math], [math]\displaystyle{ q_{n} }[/math] increases as [math]\displaystyle{ H }[/math] becomes small.
\begin{figure}[tbh] \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]\displaystyle{ \diamondsuit }[/math] indicate for [math]\displaystyle{ \omega=0.1,1.0,10 }[/math] respectively. The water depth is [math]\displaystyle{ H=20\pi }[/math] (deep water).}
(fig:rsd)
\end{center} \end{figure}
The surface displacement for line forcing is given by the sum
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]\displaystyle{ r=0 }[/math]. From the equation of motion ((2-41)) and [math]\displaystyle{ \mathbf{v}=\nabla\phi }[/math], the energy crossing a surface [math]\displaystyle{ S }[/math] in a time period [math]\displaystyle{ \left[ t,t+T\right] }[/math] is
For time-harmonic waves, i.e., [math]\displaystyle{ \phi=Re \left[ \phi\exp\left( i\omega t\right) \right] \lt math\gt , and }[/math]T=2\pi/\omega</math>, the above equation becomes
We notice that if [math]\displaystyle{ S }[/math] is the closed surface of domain [math]\displaystyle{ \mathcal{V} }[/math], then from Green's theorem, we have
which states that there is no net energy propagation to infinity, i.e., the law of energy conservation. Hence, only the travelling mode [math]\displaystyle{ H_{0}^{\left( 1\right) }\left( q_{=T= }r\right) \lt math\gt 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]\displaystyle{ 1/\sqrt{r} }[/math] and the integral on the cylindrical surface gives [math]\displaystyle{ 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]\displaystyle{ \mathrm{i}R_{n} }[/math] in Eqn.~((3-18)) and Eqn.~((3-19)) is real for each [math]\displaystyle{ 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
for the point load and
for the line load.
Note that the zeros of [math]\displaystyle{ d\left( \gamma\right) \gamma\tanh\gamma H }[/math]\ are the same as those of the dispersion equation. Thus the singularities of [math]\displaystyle{ \hat {\phi}\lt math\gt and }[/math]\hat{w}[math]\displaystyle{ 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]\displaystyle{ w_{=P= } }[/math] and [math]\displaystyle{ 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]\displaystyle{ q_{=T= } }[/math], the travelling mode carries the wave energy outwards, and the modes of [math]\displaystyle{ q_{=D= } }[/math] and [math]\displaystyle{ \mathrm{i}q_{n} }[/math], [math]\displaystyle{ n=1,2,... }[/math], the damped travelling and evanescent modes, are exponentially decaying.
====Summary of the analytic structure of [math]\displaystyle{ \hat{w====\left( \gamma\right) }[/math]}
Here we summarize and clarify the analytic properties of the complex valued function [math]\displaystyle{ \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]\displaystyle{ \hat{w} }[/math] given by Eqn.~((ap-20)) has an extra term [math]\displaystyle{ \hat{w}\left( 0\right) =0 }[/math] and is unique. Hence, even though a function such as [math]\displaystyle{ \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]\displaystyle{ \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]\displaystyle{ d\left( \gamma\right) =0 }[/math]. In addition to the residues and poles, we require that [math]\displaystyle{ \hat{w}\left( 0\right) =0 }[/math], to uniquely reconstruct [math]\displaystyle{ \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]\displaystyle{ \hat{w}\left( \gamma\right) }[/math] for [math]\displaystyle{ \gamma\in C }[/math] can be reconstructed from [math]\displaystyle{ q\in K^{} }[/math] and [math]\displaystyle{ R\left( q\right) }[/math] (either Eqn.~((3-16)) or Eqn.~((3-17))) as
However, we find the following
We note that near [math]\displaystyle{ \gamma=q\in K^{} }[/math] the left and the right hand sides share the same analytic property, i.e.,
where [math]\displaystyle{ C_{\varepsilon}\left( q\right) }[/math] is a circular contour of small radius [math]\displaystyle{ \varepsilon }[/math] around a pole [math]\displaystyle{ q }[/math]. The reason for ((3-80)) is because the function [math]\displaystyle{ 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]\displaystyle{ \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]\displaystyle{ \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]\displaystyle{ \mathrm{i}q_{n} }[/math] become equally spaced [math]\displaystyle{ n\pi/H }[/math] as the depth [math]\displaystyle{ 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]\displaystyle{ H }[/math] tends to infinity. This inspires us to find an integral expression of the summation when [math]\displaystyle{ 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]\displaystyle{ H\rightarrow\infty }[/math]. In the limit [math]\displaystyle{ \left\{ \mathrm{i}q_{n}\right\} }[/math] form a continuum with equal density over the imaginary axis. If the residue [math]\displaystyle{ R\left( \mathrm{i}q_{n}\right) }[/math] decreases proportional to [math]\displaystyle{ 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]\displaystyle{ H\approx\infty }[/math], the solution will be shown to be qualitatively unaffected by setting [math]\displaystyle{ m=0 }[/math] for typical values of [math]\displaystyle{ m }[/math] for sea ice, i.e., smaller than [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ H }[/math] becomes large. The residue at a root of the dispersion equation in Eqn.~((3-17)) tends to
as [math]\displaystyle{ H\rightarrow\infty }[/math], where
The sum over the imaginary roots in the response to point loading in Eqn.~((3-18)) is then
Similarly, the sum over imaginary roots for line load in Eqn.~((3-19)) is
The integrals in Eqn.~((3-62)) and Eqn.~((3-63)) can be evaluated by first writing [math]\displaystyle{ Q }[/math] as a sum over simple poles, as we did for the inverse Fourier transform. We notice that [math]\displaystyle{ Q\left( q\right) =\left( v\left( q\right) -v\left( -q\right) \right) /2 }[/math] where
Then [math]\displaystyle{ \omega\neq0 }[/math], [math]\displaystyle{ v\left( q\right) }[/math] and [math]\displaystyle{ v\left( -q\right) }[/math] have no poles in common and hence poles of [math]\displaystyle{ v }[/math] are also poles of [math]\displaystyle{ Q }[/math]. It follows that all poles of [math]\displaystyle{ Q }[/math] are plus and minus the poles of [math]\displaystyle{ v }[/math]. Since [math]\displaystyle{ v }[/math] and [math]\displaystyle{ w }[/math] (in the limit [math]\displaystyle{ H\rightarrow\infty }[/math]) coincide for arguments with positive real part, the poles of [math]\displaystyle{ \hat{w} }[/math] with positive real part, i.e. [math]\displaystyle{ q_{=T= },q_{=D= }\lt math\gt , and }[/math]q_{=D= }^{}[math]\displaystyle{ , are also roots of }[/math]v</math> and hence [math]\displaystyle{ Q }[/math]. There are two further roots of [math]\displaystyle{ v }[/math] with negative real part, which we denote [math]\displaystyle{ q_{=E= } }[/math] and [math]\displaystyle{ q_{=E= }^{} }[/math], with [math]\displaystyle{ q_{=E= } }[/math] chosen to have positive imaginary part. It can be shown that [math]\displaystyle{ Im \left( q_{=E= }\right) \gt 0\lt math\gt for any finite mass density }[/math]m</math>, and hence [math]\displaystyle{ 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]\displaystyle{ }[/math],[math]\displaystyle{ \triangle }[/math] and [math]\displaystyle{ \diamondsuit }[/math] indicate [math]\displaystyle{ q_{=T= } }[/math], [math]\displaystyle{ q_{=D= } }[/math] ([math]\displaystyle{ q_{=D= }^{} }[/math]) and [math]\displaystyle{ q_{=E= } }[/math] ([math]\displaystyle{ q_{=E= }^{} }[/math]).}
(fig:3-25)
\end{center} \end{figure}
Let [math]\displaystyle{ 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]\displaystyle{ v }[/math]. The residue of [math]\displaystyle{ v }[/math] at a pole [math]\displaystyle{ q\in K_{v} }[/math] is
The residues at [math]\displaystyle{ q_{=T= } }[/math], [math]\displaystyle{ q_{=D= } }[/math] and [math]\displaystyle{ q_{=D= } ^{*} }[/math], are the same as defined in section 3.2.4. We respectively denote these [math]\displaystyle{ R_{=T= } }[/math], [math]\displaystyle{ R_{=D= } }[/math] and [math]\displaystyle{ R_{=D= }^{} }[/math], as before, respectively. Write [math]\displaystyle{ R_{=E= }=R_{v}\left( q_{=E= }\right) }[/math] and hence [math]\displaystyle{ R_{=E= }^{}=R_{v}\left( q_{=E= }^{}\right) }[/math]. The residues of [math]\displaystyle{ Q }[/math] are [math]\displaystyle{ 1/2 }[/math] the residue of [math]\displaystyle{ v }[/math].
The integral in Eqn.~((3-62)) may be evaluated using the fractional function expansion
and the integral (Abramowitz and Stegun Stegun formula 11.4.47 with [math]\displaystyle{ \nu=0 }[/math], [math]\displaystyle{ t=k }[/math], [math]\displaystyle{ r=a }[/math]),
which holds for [math]\displaystyle{ Re q\gt 0 }[/math] since we always take [math]\displaystyle{ r\gt 0 }[/math]. Here [math]\displaystyle{ \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]\displaystyle{ 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]\displaystyle{ v }[/math] as,
The pole [math]\displaystyle{ -q_{=E= } }[/math] of [math]\displaystyle{ Q }[/math] has been used since [math]\displaystyle{ Re \left( -q_{=E= }\right) \gt 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
and the integral (Abramowitz and Stegun Stegun formula 5.2.13 with [math]\displaystyle{ t=k/q }[/math] and hence [math]\displaystyle{ z=qr }[/math] and formula 5.2.7)
holding for [math]\displaystyle{ Re \left( q\right) \gt 0 }[/math] since [math]\displaystyle{ r }[/math] is positive real. Here [math]\displaystyle{ Ci }[/math] and [math]\displaystyle{ 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]\displaystyle{ 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,
where, again, the pole [math]\displaystyle{ -q_{=E= } }[/math] has been used.
We emphasize here that the denominator of the function [math]\displaystyle{ 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
where [math]\displaystyle{ sgn \left( Re \gamma\right) }[/math] is sign-function of the real part of [math]\displaystyle{ \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]\displaystyle{ \gamma sgn \left( Re \gamma\right) }[/math] is often defined using the limit of [math]\displaystyle{ \sqrt{\gamma^{2}\pm i\varepsilon^{2}} }[/math] as [math]\displaystyle{ \varepsilon \rightarrow0 }[/math], which is defined on a two-sheeted Riemann surface and denoted by [math]\displaystyle{ \left| \gamma\right| _{\pm} }[/math] because it is equal to [math]\displaystyle{ \left| \gamma\right| }[/math] on the real axis.
Computation of the solutions
Static load
We can find the displacement [math]\displaystyle{ w_{=P= } }[/math] and [math]\displaystyle{ w_{=L= } }[/math]\ for static loading by setting [math]\displaystyle{ \omega=0 }[/math] in Eqn.~((3-9)), to give
and following the same procedure shown in section 3.2. The Fourier transform of the displacement is
which has four poles, [math]\displaystyle{ \pm e^{\mathrm{i}\pi/4} }[/math], [math]\displaystyle{ \pm e^{\mathrm{i}3\pi/4}\lt math\gt , 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]\displaystyle{ . The residue of }[/math]\hat{w}\left( \gamma\right) [math]\displaystyle{ at a pole }[/math]q</math> is
Hence, the solution is, from Eqn.~((3-14))
where [math]\displaystyle{ 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)
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
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]\displaystyle{ \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]\displaystyle{ 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]\displaystyle{ K_{0}\left( \zeta\right) }[/math] has the polynomial form
where [math]\displaystyle{ I_{0} }[/math] and [math]\displaystyle{ \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]\displaystyle{ \left| \zeta\right| }[/math], [math]\displaystyle{ K_{0}\left( \zeta\right) \approx-\log\zeta+c\lt math\gt where }[/math]c=\log2+\psi\left( 1\right) </math> since [math]\displaystyle{ I_{0}\left( 0\right) =1 }[/math]. Hence, as [math]\displaystyle{ r\rightarrow0 }[/math], the infinite series of [math]\displaystyle{ w_{=P= } }[/math] becomes a series of log-function of [math]\displaystyle{ r }[/math],
We notice that the second term becomes singular as [math]\displaystyle{ r\rightarrow0 }[/math].
Consider now a contour integration of the function [math]\displaystyle{ \hat{w}\left( \gamma\right) \gamma }[/math], anti-clockwise along the contour shown in Fig.~((fig:3-15)). The arc of radius [math]\displaystyle{ 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]\displaystyle{ \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]\displaystyle{ A\rightarrow\infty }[/math], [math]\displaystyle{ \hat{w}\left( \gamma\right) \gamma }[/math] tends to zero faster than [math]\displaystyle{ A^{-2} }[/math] on the semi-circle of radius [math]\displaystyle{ A }[/math], and the integral over the semi-circle tends to zero as [math]\displaystyle{ A\rightarrow\infty }[/math]. Hence the integral over the whole contour tends to zero as [math]\displaystyle{ A\rightarrow\infty }[/math]. Since this limit equals a constant multiplied by the sum of the residues of [math]\displaystyle{ \hat{w}\left( \gamma\right) \gamma }[/math] at the poles enclosed within the contour, the sum of residues of [math]\displaystyle{ \hat{w}\left( \gamma\right) \gamma }[/math] at poles in the upper half plane is zero, i.e.,
We immediately see that the term multiplied by [math]\displaystyle{ \left( c-\log r\right) /\pi }[/math] in Eqn.~((3-42)), which is the singular part, is zero. Thus at [math]\displaystyle{ r=0 }[/math] the complex displacement takes the finite value
Since [math]\displaystyle{ qR\left( q\right) }[/math] decreases as [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ H=20\pi }[/math], [math]\displaystyle{ 2\pi }[/math], and [math]\displaystyle{ 0.2\pi }[/math]. The depths are taken as multiples of [math]\displaystyle{ 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]\displaystyle{ H=2\pi\times10 }[/math], 2[math]\displaystyle{ \pi\times1 }[/math], 2[math]\displaystyle{ \pi\times0.1 }[/math], are shown. In all cases we take [math]\displaystyle{ m=0 }[/math].}
(fig:3-26)
\end{center} \end{figure}
Water depths greater than [math]\displaystyle{ 20\pi }[/math] give visually identical deflections, so the curve for [math]\displaystyle{ H=20\pi }[/math] may be taken as the deep-water solution. That is the solution for [math]\displaystyle{ H=2\pi }[/math] is nearly identical to the deep-water solution at all frequencies, so that [math]\displaystyle{ H=2\pi }[/math] may be considered deep for point load.
We notice by spliting Eqn.~((3-45)) into the real and imaginary parts,
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]\displaystyle{ H\rightarrow\infty }[/math] as described in the previous section. The sum over evanescent modes is
Using the expansion in Eqn.~((eqn:ikQexpn)) and the integral identity (Erdelyi bateman2 section 14.2 (24) with [math]\displaystyle{ a=-\mathrm{i}q }[/math], and [math]\displaystyle{ y=\mathrm{i}q }[/math])
for [math]\displaystyle{ Re q\gt 0 }[/math],\ we find that
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]\displaystyle{ r=0 }[/math] directly to give
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]\displaystyle{ H=20\pi }[/math], [math]\displaystyle{ 2\pi }[/math], and [math]\displaystyle{ 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]\displaystyle{ H=2\pi\times10 }[/math], [math]\displaystyle{ 2\pi\times1 }[/math], [math]\displaystyle{ 2\pi\times0.1 }[/math]. We take [math]\displaystyle{ m=0 }[/math].}
(fig:3-16)
\end{center} \end{figure}
We see little difference in Fig.~((fig:3-16)) between the graphs for [math]\displaystyle{ H=2\pi }[/math] and [math]\displaystyle{ 20\pi }[/math], thus [math]\displaystyle{ H }[/math] greater than [math]\displaystyle{ 2\pi }[/math] may be considered as deep for line forcing. However, at the frequency lower than [math]\displaystyle{ 1 }[/math] the graphs show difference between [math]\displaystyle{ H=2\pi }[/math] and [math]\displaystyle{ 20\pi }[/math], hence [math]\displaystyle{ 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]\displaystyle{ Ci \left( qr\right) }[/math] in terms of the log plus power series and using the identity [math]\displaystyle{ \sum_{q\in K_{v}}R_{v}\left( q\right) =0 }[/math] allows us to write
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]\displaystyle{ K_{0} }[/math] and Hankel function [math]\displaystyle{ H_{0}^{\left( 1\right) } }[/math] have a log-like singularity at [math]\displaystyle{ 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]\displaystyle{ r=0 }[/math] using a power series expansion of [math]\displaystyle{ I_{0} }[/math] given in Abramowitz and Stegun Stegun formula 9.6.10 with [math]\displaystyle{ \nu=0 }[/math]. We write the deflection as the sum of the deflection at [math]\displaystyle{ r=0 }[/math] and the remainder,
for all [math]\displaystyle{ r }[/math]. Note that near the origin the deflection behaves as [math]\displaystyle{ \left( \log\left( qr\right) +c\right) r^{2} }[/math] and thus the first derivative of the deflection near the origin is
Therefore, we have [math]\displaystyle{ 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]\displaystyle{ r }[/math] for each term in Eqn.~((3-65)) with [math]\displaystyle{ l\geq1 }[/math] has the [math]\displaystyle{ r }[/math]-dependence
which is non-zero at [math]\displaystyle{ r=0 }[/math] only for the [math]\displaystyle{ l=1 }[/math] term. Since [math]\displaystyle{ \sum_{q\in K^{}}R\left( q\right) q^{3}=1/2 }[/math], evaluating the integral of [math]\displaystyle{ \hat{w}\left( \gamma\right) \gamma^{3} }[/math] on the contour shown in Fig.~((fig:3-15)), the second derivative of the [math]\displaystyle{ l=1 }[/math] term in Eqn.~((3-65)) can be written
for small [math]\displaystyle{ r }[/math]. So we see that the strain has a singularity at the origin that behaves like [math]\displaystyle{ \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,
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
for large [math]\displaystyle{ \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
for large [math]\displaystyle{ r }[/math].
\begin{figure} [tbh] \begin{center} \includegraphics[height=7.0973cm,width=10.2297cm]{maxamp.eps} \caption{Curves of the coefficients of [math]\displaystyle{ w_{=P= } }[/math] (solid) and [math]\displaystyle{ 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
for large [math]\displaystyle{ \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]\displaystyle{ \omega=0.81 }[/math] and that of line forcing at [math]\displaystyle{ \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]\displaystyle{ w\left( r,\omega\right) }[/math] of point load at various non-dimensional frequencies [math]\displaystyle{ \omega=0.2,1.0 }[/math] and [math]\displaystyle{ 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]\displaystyle{ w\left( r,\omega\right) }[/math] for the line load problem at various non-dimensional frequencies [math]\displaystyle{ \omega=0.2,1.0 }[/math] and [math]\displaystyle{ 5.0 }[/math]. The displacement at [math]\displaystyle{ \omega=0.2 }[/math] for [math]\displaystyle{ H=2\pi }[/math] is plotted, since at a higher frequency, the response is indistinguishable between [math]\displaystyle{ H=2\pi }[/math] and [math]\displaystyle{ \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]\displaystyle{ r }[/math]. The direct formula given by Eqn.~((3-72)) is also stable enough for small [math]\displaystyle{ 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]\displaystyle{ \omega=0.2,1.0,5.0 }[/math] for the deep-water case. At low frequency [math]\displaystyle{ \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]\displaystyle{ H=2\pi }[/math] and [math]\displaystyle{ \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]\displaystyle{ l_{=c= } }[/math] and [math]\displaystyle{ 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
and the physical solution is
where [math]\displaystyle{ \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]\displaystyle{ l_{=c= } }[/math] to be [math]\displaystyle{ \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]\displaystyle{ 1 }[/math] Newton
Note that [math]\displaystyle{ \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
Therefore, for the line load the physical displacement for the line forcing, we have
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]\displaystyle{ 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]\displaystyle{ l_{=c= } }[/math]\ and normalizing scheme could be applied to the dynamic ice sheet problem. Fortunately, we find that the same characteristic length [math]\displaystyle{ l_{=c= } }[/math]\ and time [math]\displaystyle{ t_{=c= } }[/math] may be used to scale the system of equations for the range of ice thickness [math]\displaystyle{ h }[/math] and forcing frequency [math]\displaystyle{ \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]\displaystyle{ q_{=T= } }[/math], [math]\displaystyle{ q_{=D= } }[/math] and [math]\displaystyle{ q_{=E= } }[/math] of the normalized deep water solution for different ice thickness [math]\displaystyle{ 0.1 }[/math], [math]\displaystyle{ 1.0 }[/math] and [math]\displaystyle{ 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]\displaystyle{ \omega }[/math] plotted in a loglog scale. The thickness of the ice sheet is taken to be [math]\displaystyle{ 0.1 }[/math], [math]\displaystyle{ 1 }[/math] and [math]\displaystyle{ 10 }[/math] metres.}
(fig:3-5)
\end{center} \end{figure}
By closer observation of Fig.~((fig:3-5)), we find that up to [math]\displaystyle{ \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]\displaystyle{ 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]\displaystyle{ \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]\displaystyle{ H=0.2\pi }[/math] and [math]\displaystyle{ \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]\displaystyle{ H=0.2\pi }[/math] and (b) is when [math]\displaystyle{ H=\pi }[/math].}
(fig:3-17)
\end{center} \end{figure}
So far we have computed the deflection of the ice sheet with [math]\displaystyle{ m=0 }[/math]. This is justified because for typical ice sheets, it is reasonable to say that [math]\displaystyle{ m }[/math] is smaller than [math]\displaystyle{ 0.1 }[/math] assuming that ice thickness is in the order of metres and effective Young's modulus is in the order of [math]\displaystyle{ 10^{9} }[/math] Pascals, and thus the characteristic length [math]\displaystyle{ l_{=c= } }[/math] is typically larger than [math]\displaystyle{ 10h^{3/4} }[/math]. Furthermore [math]\displaystyle{ m }[/math] that is smaller than [math]\displaystyle{ 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]\displaystyle{ m }[/math] is set zero and [math]\displaystyle{ 0.1 }[/math]. The radial frequency is ranging from [math]\displaystyle{ 0.1 }[/math] to [math]\displaystyle{ 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]\displaystyle{ \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]\displaystyle{ d\left( \gamma,\omega\right) }[/math], since it is determined only by the plate equation.
We consider a smoothly shaped plate denoted by [math]\displaystyle{ \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]\displaystyle{ \mathcal{V} }[/math]. Then the Fourier transform of the plate equation must integrated over the finite domain [math]\displaystyle{ \Omega }[/math], which consists of the terms to be determined by the values of [math]\displaystyle{ w }[/math] and its derivatives at the boundary of the plate [math]\displaystyle{ \partial\Omega }[/math].
Then, using the inverse Fourier transform to calculate the displacement,
Note that [math]\displaystyle{ w\left( \bar{x},\bar{y}\right) }[/math]\ calculated using Eqn.~((3-43)) is zero outside the ice sheet [math]\displaystyle{ \Omega }[/math].
Using the Fourier transform defined above and Green's theorem, we find that the Fourier transform of the plate equation becomes
where [math]\displaystyle{ b_{1} }[/math] is function of [math]\displaystyle{ \left( \bar{\alpha},\bar{k}\right) }[/math]. [math]\displaystyle{ b_{1} }[/math] is determined by the boundary values of [math]\displaystyle{ w }[/math] and its derivatives, which results Green's theorem. The Fourier transform of Laplace's equation in [math]\displaystyle{ \Omega }[/math] for [math]\displaystyle{ -\bar{H}\lt \bar{z}\lt 0 }[/math]\ becomes
where [math]\displaystyle{ 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]\displaystyle{ b_{2} }[/math]. Thus, we find the relation between [math]\displaystyle{ \hat{\bar{\phi}}_{\Omega} }[/math] and its [math]\displaystyle{ \bar{z} }[/math] derivative similar to that in the previous sections,
where [math]\displaystyle{ b_{3} }[/math] depends on the functions [math]\displaystyle{ b_{1} }[/math] and [math]\displaystyle{ b_{2} }[/math] in Eqn.~((3-34)) and Eqn.~((3-33)). Therefore the Fourier transform [math]\displaystyle{ \hat{\bar{w}}\left( \bar{\alpha},\bar{k}\right) }[/math] has the same denominator, the dispersion function
as the infinite plate solution,
The inverse Fourier transform of [math]\displaystyle{ \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]\displaystyle{ b\left( \bar{\alpha},\bar{k}\right) }[/math]. If we assume that the numerator [math]\displaystyle{ b\left( \bar{\alpha},\bar{k}\right) }[/math] has no singularities, in other words [math]\displaystyle{ \bar{w} }[/math] and its derivatives on the edge of the plate are well-behaved functions, then the inverse Fourier transform of [math]\displaystyle{ \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]\displaystyle{ 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]\displaystyle{ b }[/math] will depend only on [math]\displaystyle{ \bar{\gamma} }[/math], the amplitude of [math]\displaystyle{ \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]\displaystyle{ b\equiv1 }[/math]. We will see an example of [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ H=0.2\pi }[/math], and when the water is deep, i.e., [math]\displaystyle{ H }[/math] larger than [math]\displaystyle{ 2\pi }[/math]. We have seen that water depth larger than [math]\displaystyle{ 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]\displaystyle{ -w_{rr}\left( r,\omega\right) }[/math], plotted as a function of non-dimensional radian frequency [math]\displaystyle{ \omega }[/math] and distance [math]\displaystyle{ r }[/math] from the point of forcing in the case when the water depth is [math]\displaystyle{ H=0.2\pi }[/math]. Then, the strain per Newton denoted by [math]\displaystyle{ \bar{S}\left( \bar{r}\right) }[/math] is calculated by
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]\displaystyle{ h/\left( 2\rho gl_{=c= }^{=4= }\right) \lt math\gt , 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]\displaystyle{ H=0.2\pi }[/math], as a function of non-dimensional radial forcing fequency (log-scale) and non-dimensional distance, [math]\displaystyle{ 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]\displaystyle{ H=0.2\pi }[/math], as a function of non-dimensional radial forcing frequency (log-scale) and non-dimensional distance, [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ r\lt 1 }[/math], changes rapidly with [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ r\gtrsim 1 }[/math] the phase of strain has a minimum value at [math]\displaystyle{ \omega\approx0.9 }[/math] and the position of the minimum is an insensitive function of distance. Hence, robust measurements can be made at [math]\displaystyle{ 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]\displaystyle{ \omega }[/math]. The acceleration is calculated from the formula of the displacement at the point of forcing, Eqn.~((3-45)), [math]\displaystyle{ -\omega^{2}w\left( 0\right) }[/math]. We notice that the both magnitude and phase of the acceleration are monotonically increasing functions of [math]\displaystyle{ \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]\displaystyle{ 0.9 }[/math] at which the shallow and deep water phases coincide (at about [math]\displaystyle{ \omega=0.9 }[/math]) so that water depth does not affect much. Let [math]\displaystyle{ \bar{\omega}_{0} }[/math] denote such frequency, then we find [math]\displaystyle{ \bar{\omega}_{0} }[/math] at which the measured relative phase is [math]\displaystyle{ 0.9 }[/math], i.e.,
From [math]\displaystyle{ \bar{\omega}_{0} }[/math] we are then able to find the characteristic length using
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
We use the asymptotic formula of a Hankel function,
Again we try to locate the physical radial frequency [math]\displaystyle{ \bar{\omega} }[/math] at which the maximum tilt occurs (at unit non-dimensional radial frequency [math]\displaystyle{ \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]\displaystyle{ \omega=\bar{\omega}t_{=c= } }[/math] by lifted ice-weight ([math]\displaystyle{ 10^{3} }[/math]kg) devided by [math]\displaystyle{ l_{=c= }^{3} }[/math] and square root of non-dimensional distance [math]\displaystyle{ \sqrt{r} }[/math], [math]\displaystyle{ 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]\displaystyle{ l_{=c= } }[/math] is originally based on the static solution, together with the characteristic time [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ \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]\displaystyle{ H=2\pi }[/math] and [math]\displaystyle{ \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,
which is about [math]\displaystyle{ 20 }[/math] to [math]\displaystyle{ 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]\displaystyle{ Ci }[/math] and [math]\displaystyle{ si }[/math] are computed by the built-in functions of MatLab. Struve function and [math]\displaystyle{ Ci }[/math] and [math]\displaystyle{ 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.
