Difference between revisions of "Wave Forcing of Small Bodies"

From WikiWaves
Jump to navigationJump to search
 
(12 intermediate revisions by 2 users not shown)
Line 1: Line 1:
=Introduction=
+
==Introduction==
  
 
While large bodies on the water surface will reflect and
 
While large bodies on the water surface will reflect and
Line 39: Line 39:
 
[[Shen and Zhong 2001]]).
 
[[Shen and Zhong 2001]]).
  
=The Marchenko Wave Drift Model=
+
==The Marchenko Wave Drift Model==
  
 
We present the models for wave forcing of small floating bodies
 
We present the models for wave forcing of small floating bodies
Line 53: Line 53:
 
system is
 
system is
 
<center><math>
 
<center><math>
m\frac{d\bar{V}_{\tau}}{dt}=-mg\frac{\bar{\eta}^{\prime}}{\sqrt{1+\bar{\eta
+
m\frac{\mathrm{d}\bar{V}_{\tau}}{\mathrm{d}t}=-mg\frac{\bar{\eta}^{\prime}}{\sqrt{1+\bar{\eta
}^{\prime}{}^{2}}} + F_d\,\,\,(1) (marchenko_original)
+
}^{\prime}{}^{2}}} + F_d\,\,\,(1)  
 
</math></center>
 
</math></center>
 
where <math>\bar{V}_{\tau}</math> is the tangential velocity and
 
where <math>\bar{V}_{\tau}</math> is the tangential velocity and
Line 65: Line 65:
 
Therefore
 
Therefore
 
<center><math>\begin{matrix}
 
<center><math>\begin{matrix}
\frac{d\bar{V}_{\tau}}{dt} & =\frac{d}{dt}\left(  \bar{V}_{\bar{x}}
+
\frac{\mathrm{d}\bar{V}_{\tau}}{\mathrm{d}t} &=\frac{\mathrm{d}}{\mathrm{d}t}\left(  \bar{V}_{\bar{x}}
 
\sqrt{1+\bar{\eta}^{\prime}{}^{2}}\right)  \\
 
\sqrt{1+\bar{\eta}^{\prime}{}^{2}}\right)  \\
& =\frac{d\bar{V}_{\bar{x}}}{dt}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}+\bar
+
&=\frac{\mathrm{d}\bar{V}_{\bar{x}}}{\mathrm{d}t}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}+\bar
{V}_{\bar{x}}\frac{d}{dt}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}\\
+
{V}_{\bar{x}}\frac{\mathrm{d}}{\mathrm{d}t}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}\\
& =\frac{d\bar{V}_{\bar{x}}}{dt}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}+\bar
+
&=\frac{\mathrm{d}\bar{V}_{\bar{x}}}{\mathrm{d}t}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}+\bar
 
{V}_{\bar{x}}\frac{1}{\sqrt{1+\bar{\eta}^{\prime}{}^{2}}}\bar{\eta}^{\prime
 
{V}_{\bar{x}}\frac{1}{\sqrt{1+\bar{\eta}^{\prime}{}^{2}}}\bar{\eta}^{\prime
}\bar{\eta}^{\prime\prime}\frac{d\bar{x}}{dt}.
+
}\bar{\eta}^{\prime\prime}\frac{\mathrm{d}\bar{x}}{\mathrm{d}t}.
 
\end{matrix}</math></center>
 
\end{matrix}</math></center>
 
Substituting this in equation (1) we obtain
 
Substituting this in equation (1) we obtain
 
<center><math>
 
<center><math>
m\frac{d\bar{V}_{\bar{x}}}{dt}=-m(g+\left(  \bar{V}_{\bar{x}}\right)  ^{2}
+
m\frac{\mathrm{d}\bar{V}_{\bar{x}}}{\mathrm{d}t}=-m(g+\left(  \bar{V}_{\bar{x}}\right)  ^{2}
 
\bar{\eta}^{\prime\prime})\frac{\bar{\eta}^{\prime}}{1+\bar{\eta}^{\prime}
 
\bar{\eta}^{\prime\prime})\frac{\bar{\eta}^{\prime}}{1+\bar{\eta}^{\prime}
{}^{2}} + F_d.\,\,\,(2) (marchenko_travelling)
+
{}^{2}} + F_d.\,\,\,(2)  
 
</math></center>
 
</math></center>
 
which is the Marchenko model in the <math>\bar{x}</math> direction.  
 
which is the Marchenko model in the <math>\bar{x}</math> direction.  
  
  
==Including added mass==
+
===Including added mass===
  
 
[[Rumer et. al. 1979]] included added mass and this makes the model
 
[[Rumer et. al. 1979]] included added mass and this makes the model
 
<center><math>
 
<center><math>
m\left(  1+C_{m}\right)  \frac{dV_{x}}{dt}=-m(g+\left(  V_{x}\right)  ^{2}
+
m\left(  1+C_{m}\right)  \frac{\mathrm{d}V_{x}}{\mathrm{d}t}=-m(g+\left(  V_{x}\right)  ^{2}
 
\eta^{\prime\prime})\frac{\eta^{\prime}}{1+\eta^{\prime}{}^{2}} + F_d,\,\,\,(3)
 
\eta^{\prime\prime})\frac{\eta^{\prime}}{1+\eta^{\prime}{}^{2}} + F_d,\,\,\,(3)
(rumer_correct_added_mass)
 
 
</math></center>
 
</math></center>
  
==Drag==
+
===Drag===
  
 
So far we have not considered the drag force although in practice the drag
 
So far we have not considered the drag force although in practice the drag
Line 109: Line 108:
 
co-ordiante system.
 
co-ordiante system.
  
=The Drift Velocity for a Linear Sinusoidal Wave=
+
==The Drift Velocity for a Linear Sinusoidal Wave==
  
 
If we are going to use the slope sliding models in the context of linear wave
 
If we are going to use the slope sliding models in the context of linear wave
Line 119: Line 118:
  
 
<center><math>
 
<center><math>
m\left(  1+C_{m}\right)  \frac{dV}{dt}=-mg\eta^{\prime}+\rho
+
m\left(  1+C_{m}\right)  \frac{\mathrm{d}V}{\mathrm{d}t}=-mg\eta^{\prime}+\rho
_{w}C_{w}A_{i}\,\,\bigg|V_{w}-V\bigg|\bigg(V_{w}-V\bigg)\,\,\,(4) (linear_force)
+
_{w}C_{w}A_{i}\,\,\bigg|V_{w}-V\bigg|\bigg(V_{w}-V\bigg)\,\,\,(4)  
 
</math></center>
 
</math></center>
 
where <math>V</math> is the velocity in the tangential or <math>x</math> direction. This equation is
 
where <math>V</math> is the velocity in the tangential or <math>x</math> direction. This equation is
Line 126: Line 125:
 
the wave profile is given by a single frequency linear wave
 
the wave profile is given by a single frequency linear wave
 
<center><math>
 
<center><math>
\eta(x,t)=\frac{H}{2}\sin(k(x-ct)),\,\,\,(5) (linear_wave)
+
\eta(x,t)=\frac{H}{2}\sin(k(x-ct)),\,\,\,(5)  
 
</math></center>
 
</math></center>
 
where <math>H</math> is the wave height and <math>k</math> is the wave number. The
 
where <math>H</math> is the wave height and <math>k</math> is the wave number. The
Line 134: Line 133:
 
</math></center>
 
</math></center>
 
where we have assumed that the water depth in infinite. Substituting equation
 
where we have assumed that the water depth in infinite. Substituting equation
(4)(linear_wave)) into (5)(linear_force)) gives  
+
(4) into (5) gives  
 
<center><math>
 
<center><math>
\left(  1+C_{m}\right)  m\frac{dV}{dt}\,=-m\,g\,\frac{kA}{2}\cos(kx-\omega
+
\left(  1+C_{m}\right)  m\frac{\mathrm{d}V}{\mathrm{d}t}\,=-m\,g\,\frac{kA}{2}\cos(kx-\omega
t) (finalsystemxb)
+
t)  
+ \,\rho_{w}AC_{w}|kcA\sin(k(x-ct))-V|(kcA\sin(k(x-ct))-V)
+
+ \,\rho_{w}AC_{w}|kcA\sin(k(x-ct))-V|(kcA\sin(k(x-ct))-V). \,\,\,(6)  
 
</math></center>
 
</math></center>
  
  
==Non-dimensionalisation==
+
===Non-dimensionalisation===
  
 
We will now non-dimensionalise equation ((finalsystemxb)). This will
 
We will now non-dimensionalise equation ((finalsystemxb)). This will
Line 148: Line 147:
 
reduce the number of variables. All length parameters are non-dimensionalised
 
reduce the number of variables. All length parameters are non-dimensionalised
 
by the amplitude <math>H/2</math> and all time parameters are non-dimensionalised by
 
by the amplitude <math>H/2</math> and all time parameters are non-dimensionalised by
<math>\sqrt{{H}/{2g}}</math> so that equation ((finalsystemxb)) becomes the following system
+
<math>\sqrt{{H}/{2g}}</math> so that equation (6) becomes the following system
 
of equations
 
of equations
 
<center><math>
 
<center><math>
Line 155: Line 154:
 
\sin{\theta}-\tilde{V}\right\vert \,\left(  \omega\sin{\theta}-\tilde{V}\right)  \\
 
\sin{\theta}-\tilde{V}\right\vert \,\left(  \omega\sin{\theta}-\tilde{V}\right)  \\
 
\frac{d\theta}{d\tilde{t}}\,=\,\omega^{2}\tilde{V}-\omega
 
\frac{d\theta}{d\tilde{t}}\,=\,\omega^{2}\tilde{V}-\omega
\end{matrix}\,\,\,(6)
+
\end{matrix}\,\,\,(7)
(autonomous)
 
 
</math></center>
 
</math></center>
 
where the variables are now non-dimensional and where
 
where the variables are now non-dimensional and where
 
<center><math>
 
<center><math>
 
\sigma=\frac{\rho_{w}C_{w}H}{2m},\quad\tau=\frac{1}{1+C_{m}
 
\sigma=\frac{\rho_{w}C_{w}H}{2m},\quad\tau=\frac{1}{1+C_{m}
},\quad\omega=\sqrt{\frac{kH}{2},\quad}=and= \quad\theta=k\left(
+
},\quad\omega=\sqrt{\frac{kH}{2},\quad} \mathrm{and} \quad\theta=k\left(
 
x-ct\right)  .
 
x-ct\right)  .
 
</math></center>
 
</math></center>
Line 176: Line 174:
 
\end{matrix}
 
\end{matrix}
 
</math></center>
 
</math></center>
Equation (6)(autonomous) is an automous system of equations which is
+
Equation (7) is an automous system of equations which is
 
periodic in <math>\theta</math> with a period of <math>2\pi</math>. The systems are therefore
 
periodic in <math>\theta</math> with a period of <math>2\pi</math>. The systems are therefore
 
defined on a cylinder. It depends on three non-dimensional parameters, of
 
defined on a cylinder. It depends on three non-dimensional parameters, of
Line 187: Line 185:
 
It should also be noted that all the constants must be positive.
 
It should also be noted that all the constants must be positive.
  
==The behaviour of the solutions==
+
===The behaviour of the solutions===
  
 
Analytic solutions for system (4) cannot be obtained easily.
 
Analytic solutions for system (4) cannot be obtained easily.
Line 195: Line 193:
 
about the behaviour of the equations using ideas from dynamical systems theory.
 
about the behaviour of the equations using ideas from dynamical systems theory.
  
If the velocity of the body, <math>V,</math> is large, <math>dV/dt</math> can therefore be
+
If the velocity of the body, <math>V,</math> is large, <math>\mathrm{d}V/\mathrm{d}t</math> can therefore be
 
approximated by
 
approximated by
 
<center><math>
 
<center><math>
\frac{dV}{dt}\approx-\sigma\tau\,|V|V,
+
\frac{\mathrm{d}V}{\mathrm{d}t}\approx-\sigma\tau\,|V|V,
 
</math></center>
 
</math></center>
 
which means that if <math>V</math> is large and positive the <math>V</math> will decrease, if <math>V</math> is
 
which means that if <math>V</math> is large and positive the <math>V</math> will decrease, if <math>V</math> is
Line 222: Line 220:
 
<center><math>\begin{matrix}
 
<center><math>\begin{matrix}
 
\iint_{S}\left(  \frac{\partial P}{\partial V}+\frac{\partial Q}
 
\iint_{S}\left(  \frac{\partial P}{\partial V}+\frac{\partial Q}
{\partial\theta}\right)  \ dVd\theta &  =\oint_{\Gamma}\left(  Pd\theta
+
{\partial\theta}\right)  \ \mathrm{d}V\mathrm{d}\theta &  =\oint_{\Gamma}\left(  P\mathrm{d}\theta
-QdV\right)  \\
+
-Q\mathrm{d}V\right)  \\
&  =\int_{0}^{T}\left(  P\frac{d\theta}{dt}-Q\frac{dV}{dt}\right)  dt\\
+
&  =\int_{0}^{T}\left(  P\frac{\mathrm{d}\theta}{\mathrm{d}t}-Q\frac{\mathrm{d}V}{\mathrm{d}t}\right)  \mathrm{d}t\\
&  =\int_{0}^{T}\left(  PQ-QP\right)  dt=0,
+
&  =\int_{0}^{T}\left(  PQ-QP\right)  \mathrm{d}t=0,
 
\end{matrix}</math></center>
 
\end{matrix}</math></center>
 
which is known as Bendixon's criterion. However,
 
which is known as Bendixon's criterion. However,
Line 235: Line 233:
 
<center><math>
 
<center><math>
 
\iint_{S}\frac{\partial P}{\partial V}+\frac{\partial Q}{\partial\theta
 
\iint_{S}\frac{\partial P}{\partial V}+\frac{\partial Q}{\partial\theta
}\ dVd\theta>0
+
}\ \mathrm{d}V\mathrm{d}\theta>0
 
</math></center>
 
</math></center>
 
and  a closed limit cycle cannot exist. Futhermore, we can use a similar
 
and  a closed limit cycle cannot exist. Futhermore, we can use a similar
 
arguement to show that there can be at most on limit cycle which encircles the
 
arguement to show that there can be at most on limit cycle which encircles the
cylinder. This means that the solution to equation ((autonomous)) and
+
cylinder. This means that the solution to equation (7) and
hence to equation ((finalsystemxb)) much tend to either an equilibrium
+
hence to equation (6) much tend to either an equilibrium
 
point or a encircling limit cycle. The equilibrium points are characterised by
 
point or a encircling limit cycle. The equilibrium points are characterised by
 
a solution which "surfs" the wave, i.e. it stays at a fixed phase of the wave
 
a solution which "surfs" the wave, i.e. it stays at a fixed phase of the wave
Line 247: Line 245:
 
is false (as long as the surface friction is not assumed to be zero).
 
is false (as long as the surface friction is not assumed to be zero).
  
==Equilibrium points==
+
===Equilibrium points===
  
 
We will now investigate the existence of the equilibrium points. For this
 
We will now investigate the existence of the equilibrium points. For this
Line 255: Line 253:
 
valid. At an equilibrium point <math>(\theta_{j},V_{j})</math> the conditions
 
valid. At an equilibrium point <math>(\theta_{j},V_{j})</math> the conditions
 
<center><math>\begin{matrix}
 
<center><math>\begin{matrix}
\frac{dV}{dt} &  =0,\\
+
\frac{\mathrm{d}V}{\mathrm{d}t} &  =0,\\
\frac{d\theta}{dt} &  =0,
+
\frac{\mathrm{d}\theta}{\mathrm{d}t} &  =0,
 
\end{matrix}</math></center>
 
\end{matrix}</math></center>
 
are satisfied. The second equation implies <math>V_{j}=1/\omega,</math> (which means that
 
are satisfied. The second equation implies <math>V_{j}=1/\omega,</math> (which means that
Line 263: Line 261:
 
<center><math>
 
<center><math>
 
-\omega^{2}\,\cos{\theta}+\sigma\left\vert \omega\sin{\theta}-1/\omega
 
-\omega^{2}\,\cos{\theta}+\sigma\left\vert \omega\sin{\theta}-1/\omega
\right\vert \,\left(  \omega\sin{\theta}-1/\omega\right)  =0
+
\right\vert \,\left(  \omega\sin{\theta}-1/\omega\right)  =0\,\,\,(8)
(equipoint)
 
 
</math></center>
 
</math></center>
 
If <math>0\leq\theta\leq{\pi}/{2}</math> or <math>{3\pi}/{2}\leq\theta\leq2\pi</math>, this equation
 
If <math>0\leq\theta\leq{\pi}/{2}</math> or <math>{3\pi}/{2}\leq\theta\leq2\pi</math>, this equation
Line 287: Line 284:
 
<math>0\leq\theta\leq{\pi}/{2}</math> or <math>{3\pi}/{2}\leq\theta\leq2\pi</math> it follows that
 
<math>0\leq\theta\leq{\pi}/{2}</math> or <math>{3\pi}/{2}\leq\theta\leq2\pi</math> it follows that
 
<center><math>
 
<center><math>
 
+
\frac{\partial P}{\partial\theta}(\theta_{1},0) > 0,\quad \mathrm{and}\quad
\frac{\partial P}{\partial\theta}(\theta_{1},0)\  &  >\ 0,\\
+
\frac{\partial P}{\partial\theta}(\theta_{2},0)  < 0.\,\,\,(9)
\frac{\partial P}{\partial\theta}(\theta_{2},0)\  & <\ 0.
 
 
 
(equiclassi)
 
 
</math></center>
 
</math></center>
 
We can determine the type of the equilibrium points by considering the
 
We can determine the type of the equilibrium points by considering the
Line 298: Line 292:
 
\left(
 
\left(
 
\begin{matrix}
 
\begin{matrix}
[c]{cc}
 
 
-2\sigma\left\vert \omega\sin\omega-1/\omega\right\vert  & \frac{\partial
 
-2\sigma\left\vert \omega\sin\omega-1/\omega\right\vert  & \frac{\partial
 
P}{\partial\theta}\\
 
P}{\partial\theta}\\
Line 312: Line 305:
 
^{2}\theta}}}.
 
^{2}\theta}}}.
 
</math></center>
 
</math></center>
It follows from ((equiclassi)) that at <math>\theta_{1}</math> there is one
+
It follows from (9) that at <math>\theta_{1}</math> there is one
 
eigenvalue with positive real part and one with negative real part. At
 
eigenvalue with positive real part and one with negative real part. At
 
<math>\theta_{2}</math> both eigenvalues have a negative real part. Applying the
 
<math>\theta_{2}</math> both eigenvalues have a negative real part. Applying the
equilibrium point classification theorems ([[perko91]]), <math>\theta_{1}</math> is a
+
equilibrium point classification theorems, <math>\theta_{1}</math> is a
 
saddle point and <math>\theta_{2}</math> is an attracting node or a spiral.  
 
saddle point and <math>\theta_{2}</math> is an attracting node or a spiral.  
  
In summary we have shown that all solutions to equation ((autonomous)) and
+
In summary we have shown that all solutions to equation (7) and
hence to equation ((finalsystemxb)) must tend to either an equilibrium
+
hence to equation (6) must tend to either an equilibrium
 
point (at which the velocity of the body is given by the wave phase speed) or
 
point (at which the velocity of the body is given by the wave phase speed) or
 
to an encircling limit cycle. There can exist at most one encircling limit
 
to an encircling limit cycle. There can exist at most one encircling limit
Line 325: Line 318:
 
small and they come in an attracting node and saddle pair (a saddle-node
 
small and they come in an attracting node and saddle pair (a saddle-node
 
bifurcation).
 
bifurcation).
 
=Drift Vevolcity: Numerical solution=
 
 
We now investigate the numerical solution of equation ((autonomous)). It
 
turns out, that although it is easy to solve this equation numerically there
 
are some sublte problems associated with determining the drift velocity.
 
 
Numerical solutions are presented for two settings, one in which all solutions
 
tend to a limit cycle and one in which there also exist equilibrium points. In
 
addition to this, a numerical procedure to calculate the steady state solution
 
is described. Using this method, the long term drift velocity is calculated
 
and its dependence on the different parameters is investigated.
 
 
==Phase plots==
 
 
The autonomous system of differential equations (autonomous) can be
 
solved numerically using standard ode-solving methods, for example the
 
Runge-Kutta 5(4) method.
 
 
The system depends on three non-dimensionalised parameters <math>\omega</math> <math>\tau</math>,
 
and <math>
 
\sigma<math>, which have to be determined. </math>\omega</math> can be calculated directly from
 
the wavelength and the amplitude of the wave by <math>\omega= \sqrt{2 \pi A/L }</math>.
 
To determine <math>\sigma</math>, the relative density of the ice to that of the water,
 
the non-dimensionalised height <math>h</math> of the ice floe (which is the fraction of
 
height of the ice floe and amplitude of the wave) and the drag coefficient
 
have to be known. We will begin with the values used in [[Shen and Zhong 2001]].
 
Shen \& Zhong assumed that the added mass was <math>C_m=0.08</math> so
 
<math>\tau=1/(1.08)=0.93</math>.
 
They assumed that the floating body was of constant cross section so
 
that the mass was <math>m=\rho_i A D</math> where <math>\rho_i</math> is the density
 
and <math>D</math> is the body thickness. This means that we can express <math>\sigma</math>
 
as
 
<center><math>\sigma =\frac{\rho_w C_w H}{2\rho_i D}.</math></center>
 
Shen \& Zhong further assumed that <math>H/\lambda = 0.02</math>, that
 
<math>C_w\lambda/D = 0.1</math> where <math>\lambda</math> is the wavelength and that
 
<math>\rho_w/\rho_i = 1/0.9</math> This gives
 
us a value for <math>\sigma</math> of <math>\sigma = 0.02/1.8=0.011</math>. The value
 
of <math>\omega</math> is found from <math>\omega = \sqrt{\pi H/\lambda}</math> so
 
that <math>\omega = \sqrt{2\pi\times 0.01} = 0.25</math>. We can see that
 
<math>\omega^2/(\omega-1/\omega)^2 = 0.0045</math> which is less than
 
<math>\sigma</math> so we do not have any
 
equilibrium points.
 
 
Rike - I want the plot with the values above. Plus the
 
same values except now <math>\sigma = 0.002/1.8</math> which
 
should now have equilibrium points. If it does
 
not then email me and we will find closely related
 
values which do.
 
 
Figures (mymodelvec04) and (mymodelvec08) show phase plots for
 
different settings for the wavelength and the amplitude of the wave. The plots
 
are projections of a cylinder, which means that there is a mapping from the
 
right edge of the picture to the left edge.
 
 
\begin{figure}[p]
 
\begin{center}
 
\includegraphics[width=9cm]{mymodelvec04}
 
\end{center}
 
\par
 
 
 
(mymodelvec04)\end{figure}
 
 
\begin{figure}[p]
 
\begin{center}
 
\includegraphics[width=9cm]{mymodelvec08}
 
\end{center}
 
\par
 
 
 
(mymodelvec08)\end{figure}
 
 
Figure (mymodelvec04) illustrates the situation where the wavelength is 40
 
times bigger than the amplitude of the wave -- the wave is very flat. In this
 
case <math>\omega</math> is <math>0.4</math>. All solutions tend towards one attractive limit cycle
 
and no equilibrium points exist.\newline In figure (mymodelvec08) it is
 
assumed that the wavelength is approximately 10 times bigger than the wave
 
amplitude, so that <math>\omega</math> is given by 0.8. In this case, two equilibrium
 
points, <math>(\theta_{1}, {1}/{\omega}) </math> and <math>(\theta_{2}, {1}/{\omega})</math>, and
 
one limit cycle exist. Therefore the cylinder can be divided into two regions:
 
One in which all solutions tend towards the limit cycle, and another in which
 
all solutions tend to the attractive equilibrium point. The border between
 
these regions is given by the two solutions which tend to the saddle point and
 
the solution at the left side of the picture (this is the dashed solution
 
which tends towards the saddle point in the next period).
 
 
The solutions which come out of the saddle point can be found by starting
 
close to the saddle point and solving system ((autonomous)). The others
 
are also found by starting close to the saddle point, but solving system
 
((autonomous)) where the right hand sides are multiplied by <math>-1 </math>.
 
<math>\theta_{1}</math> and <math>\theta_{2}</math> can be found by solving equation
 
((equipoint)), for example by finding the zeros of <math>P(\theta,{1}/{\omega
 
})</math> numerically with Newton's method. The starting points can be chosen as
 
<math>x_{0}={\pi}/{2}</math> and as <math>x_{0}={3\pi}/{2}</math>. It can be seen that the two
 
equilibrium points are situated close to the wave trough and the wave crest.
 
 
==The limit cycle==
 
 
Determining the limit cycle accurately is of great
 
importance for in the calculation of the drift velocity. However it is
 
not a simply calculation because of the very slow convergence of the solution
 
to the limit cycle. For example, it is not a good idea to just let the
 
solution evolve and then calculate the drift velocity and was done by
 
[[Shen and Zhong 2001]] because the convergence is so slow that it may appear as though the
 
solution has converged when it has not.
 
 
The problem which has to be solved is a boundary value problem with periodic
 
boundary conditions. For clarity, in this section <math>V(\theta;\theta_{0},v_{0})</math>
 
denotes the velocity <math>V(\theta)</math> satisfying the initial condition
 
<math>v(\theta_{0})=v_{0}</math>. It should be noted that the dependent
 
variable <math>t</math> does not appear in the above conditions which
 
greatly complicates the problem of finding the solution.
 
 
To find the limit cycle, two problems have to be solved,
 
 
 
 
#What is the initial velocity <math>v_{0}</math> at <math>\theta_{0} = 2\pi</math> such that
 
the velocity at <math>\theta= 0</math>, <math>v(0;2\pi,v_{0}) = v_{1}</math>, is equal to the
 
initial velocity, <math>v_{0} = v_{1}</math>? This involves the multiple solving of the sub-problem:
 
 
#Given an initial condition <math>(2\pi, v_{0})</math>, what is the velocity <math>v_{1}
 
= v(0;2\pi,v_{0})</math>?
 
 
 
Since the it is not possible to solve the first problem without solving the
 
second we begin by considering the second problem. One method to solve the
 
second problem is using interpolation. The system of differential equations is
 
solved, using very small timesteps to obtain good accuracy. The solution
 
vector returned by the ode-solver is searched for the two closest values to
 
<math>\theta=0</math>: The closest negative value, <math>\theta_{-}</math>, and the closest positive
 
one, <math>\theta_{+}</math>. Using these values, a good approximation of the velocity at
 
<math>\theta=0</math> and the required time can be found by linear
 
interpolation.
 
 
The easiest way to solve the first problem is
 
to use the fact that each solution will eventually tend towards
 
the limit cycle. The starting initial condition is <math>(2\pi, v_{0})</math>. The
 
velocity <math>v_{1} = v(0;2\pi,v_{0})</math> can now be found using the method described
 
above. Since the solution is periodic, this velocity <math>v_{1}</math> is used as the
 
new initial condition: The system is solved again with initial condition
 
<math>(2\pi, v_{1})</math> and the whole method is iterated. The iteration stops when
 
<math>|v_{n} - v_{n-1}|</math> is smaller than a tolerance <math>\varepsilon</math>.
 
However this iterative  method is very slow.
 
It is much better to use a  numerical
 
algorithms for finding zeros. The fixed point <math>v_{0} = v_{1}</math> can be found by
 
finding the zero of
 
<center><math>\begin{matrix}
 
(zerosof)g(v_{0}) = v_{0} - v(0;2\pi,v_{0}) = v_{0} - v_{1},
 
\end{matrix}</math></center>
 
where <math>v_{0}</math> is the initial velocity at <math>\theta_{0}=2\pi</math> and <math>v_{1}</math> is the
 
corresponding velocity at <math>\theta=0</math>. The zero of equation ((zerosof))
 
are found here by the secant method.
 
 
==Finding the Drift Velocity from the Limit Cycle==
 
 
In theory the drift velocity can be calculated from the
 
the time <math>\tilde{t}</math> the ice floe needs to travel
 
from one wavecrest to the next wavecrest using the formula
 
<center><math>\begin{matrix}
 
v_{d} = c - \frac{L}{\tilde{t}}.
 
\end{matrix}</math></center>
 
However, a small change in the value of <math>\tilde{t} </math> in this
 
calculation results in a great variation in <math>v_{d}</math>. For this reason
 
we use the an alternative method in which
 
the velocity of the ice floe with respect to time is
 
integrated over one period. This result is divideed by <math>\tilde{t}</math>
 
to obtain the drift velocity.
 
 
You need to redraw these figures using the non-dimensional variables
 
 
\begin{figure}[h]
 
\begin{center}
 
\includegraphics[width=0.9\textwidth]{mymodeldriftvsdrag}
 
\end{center}
 
\caption{Drift velocity for different drag coefficients and different
 
amplitudes}
 
(driftversusdrag)
 
\end{figure}
 

Latest revision as of 12:11, 11 December 2009

Introduction

While large bodies on the water surface will reflect and scatter waves, if the wavelength is much longer than the dimension of the body, the wavefield will be little modified.. In this case the wave diffraction will be negligible and the object will be passively driven by the waves. The study of this passive drift of floating bodies is important for predicting the drift of bouyant debris. While there is a wide range of application, the geophysical and offshore engineering problem of the wave drift of small icefloes and iceburg debris has been the motivation of much of the research in the wave induced drift of small floating bodies.

The first model for the wave induced drift of small floating bodies was developed by Rumer et. al. 1979 to model the drift of ice floes in the great lakes. This model was based on decomposing the wave force into two components, the first due to the drag between the body and water and the second due to the sliding effect of the body on the surface of the wave. The Rumer model was later used by Shen and Ackley 1991 to investigate the drift of pankcake ice. The Rumer model was then further investigated by Shen and Zhong 2001 where consideration was given to the effect of reflection and where analytic solutions were derived in limiting situations. Shen and Zhong 2001 also presented results showing that the wave drift is a function of the initial body position and velcocity. In all cases the underlying wave field was assumed to be the small amplitude linear wave.

Independently of the model developed by Rumer, Marchenko 1999 derived a model for the wave induced drift of small ice floes in waves. Like Rumer, Markenko decomposed the wave force into two components, one due to drag the other due to sliding. The equations of motion in Marchenko's formulation were in terms of the normal and tangential directions of the wave surface and therefore were difficult to compare to the similar equations of Rumer.

Grotmaack and Meylan 2006 compared the models of Rumer and Marchenko and establish that the models are not the same. By a third derivation method they established that the correct model is Marchenko's. They also showed that the long term drift velocity cannot be a function of either the initial position or initial velocity (contradicting the results for drift velocity presented in Shen and Zhong 2001).

The Marchenko Wave Drift Model

We present the models for wave forcing of small floating bodies which have been derived by Marchenko 1999. Marchenko decomposed the force acting of the small body into two components due to the drag force between the water and the body and the gravity force due to the body sliding down the surface of the wave. The drag force is due to the difference between the body and water velocities squared. In the Marchenko model the coordinate system travels with the wave and the velocity in the Marchenko model is in the tangential direction.

The model for the sliding force given by Marchenko in the moving co-ordinate system is

[math]\displaystyle{ m\frac{\mathrm{d}\bar{V}_{\tau}}{\mathrm{d}t}=-mg\frac{\bar{\eta}^{\prime}}{\sqrt{1+\bar{\eta }^{\prime}{}^{2}}} + F_d\,\,\,(1) }[/math]

where [math]\displaystyle{ \bar{V}_{\tau} }[/math] is the tangential velocity and [math]\displaystyle{ F_d }[/math] is the drag force (which we will introduce later). The velocity in the [math]\displaystyle{ \bar{x} }[/math] direction is given by

[math]\displaystyle{ \bar{V}_{\bar{x}}=\frac{\bar{V}_{\tau}}{\sqrt{1+\bar{\eta}^{\prime}{}^{2}}} }[/math]

Therefore

[math]\displaystyle{ \begin{matrix} \frac{\mathrm{d}\bar{V}_{\tau}}{\mathrm{d}t} &=\frac{\mathrm{d}}{\mathrm{d}t}\left( \bar{V}_{\bar{x}} \sqrt{1+\bar{\eta}^{\prime}{}^{2}}\right) \\ &=\frac{\mathrm{d}\bar{V}_{\bar{x}}}{\mathrm{d}t}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}+\bar {V}_{\bar{x}}\frac{\mathrm{d}}{\mathrm{d}t}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}\\ &=\frac{\mathrm{d}\bar{V}_{\bar{x}}}{\mathrm{d}t}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}+\bar {V}_{\bar{x}}\frac{1}{\sqrt{1+\bar{\eta}^{\prime}{}^{2}}}\bar{\eta}^{\prime }\bar{\eta}^{\prime\prime}\frac{\mathrm{d}\bar{x}}{\mathrm{d}t}. \end{matrix} }[/math]

Substituting this in equation (1) we obtain

[math]\displaystyle{ m\frac{\mathrm{d}\bar{V}_{\bar{x}}}{\mathrm{d}t}=-m(g+\left( \bar{V}_{\bar{x}}\right) ^{2} \bar{\eta}^{\prime\prime})\frac{\bar{\eta}^{\prime}}{1+\bar{\eta}^{\prime} {}^{2}} + F_d.\,\,\,(2) }[/math]

which is the Marchenko model in the [math]\displaystyle{ \bar{x} }[/math] direction.


Including added mass

Rumer et. al. 1979 included added mass and this makes the model

[math]\displaystyle{ m\left( 1+C_{m}\right) \frac{\mathrm{d}V_{x}}{\mathrm{d}t}=-m(g+\left( V_{x}\right) ^{2} \eta^{\prime\prime})\frac{\eta^{\prime}}{1+\eta^{\prime}{}^{2}} + F_d,\,\,\,(3) }[/math]

Drag

So far we have not considered the drag force although in practice the drag force is the more difficult force to model because it depends on an unknown factor which models the friction force between the body and the water. However there is a consensus that the drag force should be proprotional to the square of the velocity difference of the water particles at the water surface and the small body. The model for the drag force is therefore given by

[math]\displaystyle{ F_{d}=\rho_{w}C_{w}A\,\,\bigg|V_{w}-V\bigg|\bigg(V_{w}-V\bigg) }[/math]

where [math]\displaystyle{ \rho }[/math] is the density of the medium through which the body moves, [math]\displaystyle{ A }[/math] is the area of the moving object, [math]\displaystyle{ C_{w} }[/math] is the drag coefficient, [math]\displaystyle{ V }[/math] is the velocity given in the [math]\displaystyle{ x }[/math] or tangential directions as appropriate and [math]\displaystyle{ V_{w} }[/math] is the velocity of the water particles also given in the appropriate co-ordiante system.

The Drift Velocity for a Linear Sinusoidal Wave

If we are going to use the slope sliding models in the context of linear wave theory then it makes sense to work derive the system of equations under the same assumptions which underlie the linear wave theory. This means that the tangential and [math]\displaystyle{ x }[/math] directions should be considered as equivalent and that higher order terms should be neglected. Under these assumptions equation ((rumer_correct_added_mass)) becomes

[math]\displaystyle{ m\left( 1+C_{m}\right) \frac{\mathrm{d}V}{\mathrm{d}t}=-mg\eta^{\prime}+\rho _{w}C_{w}A_{i}\,\,\bigg|V_{w}-V\bigg|\bigg(V_{w}-V\bigg)\,\,\,(4) }[/math]

where [math]\displaystyle{ V }[/math] is the velocity in the tangential or [math]\displaystyle{ x }[/math] direction. This equation is identical to the equation which is used in Shen and Zhong 2001. We asume that the wave profile is given by a single frequency linear wave

[math]\displaystyle{ \eta(x,t)=\frac{H}{2}\sin(k(x-ct)),\,\,\,(5) }[/math]

where [math]\displaystyle{ H }[/math] is the wave height and [math]\displaystyle{ k }[/math] is the wave number. The velocity of a particle at the water surface is given by

[math]\displaystyle{ V_{w}=\frac{kcH}{2}\sin(k(x-ct)) }[/math]

where we have assumed that the water depth in infinite. Substituting equation (4) into (5) gives

[math]\displaystyle{ \left( 1+C_{m}\right) m\frac{\mathrm{d}V}{\mathrm{d}t}\,=-m\,g\,\frac{kA}{2}\cos(kx-\omega t) + \,\rho_{w}AC_{w}|kcA\sin(k(x-ct))-V|(kcA\sin(k(x-ct))-V). \,\,\,(6) }[/math]


Non-dimensionalisation

We will now non-dimensionalise equation ((finalsystemxb)). This will serve two purposes, the first to simplify the equations and the second to reduce the number of variables. All length parameters are non-dimensionalised by the amplitude [math]\displaystyle{ H/2 }[/math] and all time parameters are non-dimensionalised by [math]\displaystyle{ \sqrt{{H}/{2g}} }[/math] so that equation (6) becomes the following system of equations

[math]\displaystyle{ \begin{matrix} \frac{d\tilde{V}}{d\tilde{t}}\,=-\tau\omega^{2}\,\cos{\theta}+\sigma\tau\left\vert \omega \sin{\theta}-\tilde{V}\right\vert \,\left( \omega\sin{\theta}-\tilde{V}\right) \\ \frac{d\theta}{d\tilde{t}}\,=\,\omega^{2}\tilde{V}-\omega \end{matrix}\,\,\,(7) }[/math]

where the variables are now non-dimensional and where

[math]\displaystyle{ \sigma=\frac{\rho_{w}C_{w}H}{2m},\quad\tau=\frac{1}{1+C_{m} },\quad\omega=\sqrt{\frac{kH}{2},\quad} \mathrm{and} \quad\theta=k\left( x-ct\right) . }[/math]

We can think of [math]\displaystyle{ \sigma }[/math] as the non-dimensional drag coefficient, [math]\displaystyle{ \tau }[/math] as the corrected mass, [math]\displaystyle{ \omega }[/math] as the wave frequency and [math]\displaystyle{ \theta }[/math] as the co-ordinates moving with the wave. From now on we will drop the tilde and assume that all variables are non-dimensional. We will introduce the following notation which we will need later

[math]\displaystyle{ \begin{matrix} P(V,\theta)=-\omega^{2}\,\cos{\theta}+\sigma\left\vert \omega\sin{\theta }-V\right\vert \,\left( \omega\sin{\theta}-V\right) \\ Q\left( V,\theta\right) =\,\omega^{2}V-\omega \end{matrix} }[/math]

Equation (7) is an automous system of equations which is periodic in [math]\displaystyle{ \theta }[/math] with a period of [math]\displaystyle{ 2\pi }[/math]. The systems are therefore defined on a cylinder. It depends on three non-dimensional parameters, of which only [math]\displaystyle{ \sigma }[/math] and [math]\displaystyle{ \omega }[/math] are really important since [math]\displaystyle{ \tau }[/math] must be close to unity. We can write the constant [math]\displaystyle{ \omega }[/math] (using the assumption of infinite depth) as

[math]\displaystyle{ \omega=\sqrt{2\pi}\sqrt{\frac{A}{\lambda}}, }[/math]

It should also be noted that all the constants must be positive.

The behaviour of the solutions

Analytic solutions for system (4) cannot be obtained easily. Starting with the simplified Rumer et. al. 1979 model (without centripetal force), Shen and Zhong 2001 found approximate analytic solutions for a few special cases, for example when [math]\displaystyle{ C_{w}=0 }[/math]. We will derive here quantitative results about the behaviour of the equations using ideas from dynamical systems theory.

If the velocity of the body, [math]\displaystyle{ V, }[/math] is large, [math]\displaystyle{ \mathrm{d}V/\mathrm{d}t }[/math] can therefore be approximated by

[math]\displaystyle{ \frac{\mathrm{d}V}{\mathrm{d}t}\approx-\sigma\tau\,|V|V, }[/math]

which means that if [math]\displaystyle{ V }[/math] is large and positive the [math]\displaystyle{ V }[/math] will decrease, if [math]\displaystyle{ V }[/math] is large and negative [math]\displaystyle{ V }[/math] will increase. This means that the solution, which lives on the cylinder -[math]\displaystyle{ \pi\leq\theta\leq\pi }[/math], [math]\displaystyle{ -\infty\lt V\lt \infty, }[/math] cannot leave a bounded region of the cylinder. This means that we can use the Poincare-Bendixson Theorem, which in the case of a cylinder tells us that the solution will either tend to a limit cycle or towards an equilibrium point. It should be noted that there are two possible limit cycles on a cylinder, and ordinary limit cycle that can be shrunk continuously to a point and a limit cycle which encircles the cylinder itself. We refer to the first limit cycle as a closed limit cycle and the second limit cycle as an encircling limit cycle .

We can easily establish that there cannot be a closed limit cycle by the following argument. Suppose there exists a closed limit cycle, which we denote by [math]\displaystyle{ \Gamma, }[/math] i.e.

[math]\displaystyle{ \Gamma=\left\{ ({V(t)},{\theta(t)})\,\big|\,0\leq t\leq T\right\} }[/math]

and [math]\displaystyle{ \Gamma\left( 0\right) =\Gamma\left( T\right) }[/math] in a simply connected region of the cylinder [math]\displaystyle{ S }[/math]. From Green's theorem it follows that follows

[math]\displaystyle{ \begin{matrix} \iint_{S}\left( \frac{\partial P}{\partial V}+\frac{\partial Q} {\partial\theta}\right) \ \mathrm{d}V\mathrm{d}\theta & =\oint_{\Gamma}\left( P\mathrm{d}\theta -Q\mathrm{d}V\right) \\ & =\int_{0}^{T}\left( P\frac{\mathrm{d}\theta}{\mathrm{d}t}-Q\frac{\mathrm{d}V}{\mathrm{d}t}\right) \mathrm{d}t\\ & =\int_{0}^{T}\left( PQ-QP\right) \mathrm{d}t=0, \end{matrix} }[/math]

which is known as Bendixon's criterion. However,

[math]\displaystyle{ \frac{\partial P}{\partial V}+\frac{\partial Q}{\partial\theta}\,=\,-2\sigma \tau\,|\omega\sin{\theta}-V|, }[/math]

which is always negative except on the line [math]\displaystyle{ V=\omega\sin{\theta} }[/math]. Therefore

[math]\displaystyle{ \iint_{S}\frac{\partial P}{\partial V}+\frac{\partial Q}{\partial\theta }\ \mathrm{d}V\mathrm{d}\theta\gt 0 }[/math]

and a closed limit cycle cannot exist. Futhermore, we can use a similar arguement to show that there can be at most on limit cycle which encircles the cylinder. This means that the solution to equation (7) and hence to equation (6) much tend to either an equilibrium point or a encircling limit cycle. The equilibrium points are characterised by a solution which "surfs" the wave, i.e. it stays at a fixed phase of the wave and travels at the wave phase speed. It is clear now that the claim in Shen and Zhong 2001 that the drift velocity depends on the initial conditions is false (as long as the surface friction is not assumed to be zero).

Equilibrium points

We will now investigate the existence of the equilibrium points. For this section, we will make a futher assumption that [math]\displaystyle{ \omega\lt 1 }[/math]. This assumption will be valid for all but the steepest waves and for waves so steep that [math]\displaystyle{ \omega\geq1 }[/math] then the linear assumptions we have made will no longer be valid. At an equilibrium point [math]\displaystyle{ (\theta_{j},V_{j}) }[/math] the conditions

[math]\displaystyle{ \begin{matrix} \frac{\mathrm{d}V}{\mathrm{d}t} & =0,\\ \frac{\mathrm{d}\theta}{\mathrm{d}t} & =0, \end{matrix} }[/math]

are satisfied. The second equation implies [math]\displaystyle{ V_{j}=1/\omega, }[/math] (which means that at the equilibrium point the velocity is fixed to be the wave phase speed as expected) and if we substitute this equation into the first we obtain.

[math]\displaystyle{ -\omega^{2}\,\cos{\theta}+\sigma\left\vert \omega\sin{\theta}-1/\omega \right\vert \,\left( \omega\sin{\theta}-1/\omega\right) =0\,\,\,(8) }[/math]

If [math]\displaystyle{ 0\leq\theta\leq{\pi}/{2} }[/math] or [math]\displaystyle{ {3\pi}/{2}\leq\theta\leq2\pi }[/math], this equation cannot be satisfied since then the left hand side is always negative (remember that we have assumed that [math]\displaystyle{ \omega\lt 1 }[/math]. We can also see that a necessary condition for equilibruim points is that

[math]\displaystyle{ \sigma\lt \frac{\omega^{2}}{\left( \omega-1/\omega\right) ^{2}} }[/math]

(again using our assumption that [math]\displaystyle{ \omega\lt 1). }[/math] This makes sense, because for a body to be moving at the speed of the wave the drag force must be small compared to the sliding force. A graphical analysis shows that at most two equilibrium points can exist and that, for given [math]\displaystyle{ \sigma }[/math], [math]\displaystyle{ \omega }[/math] must be large enough to allow the existence of equilibrium points. Therefore, for a given drag there is a frequency (or wave height) below which no equlilibrium points exist. In most practical situation there are no equlilibrium points (which explains why they were not observed in Shen and Zhong 2001. We will denote the two equilibrium points by [math]\displaystyle{ \theta_{1} }[/math] and [math]\displaystyle{ \theta_{2} }[/math] (and not consider the critical case where there is only on equilibrium point) and assume that [math]\displaystyle{ \theta_{1} }[/math] is smaller than [math]\displaystyle{ \theta_{2} }[/math] so that [math]\displaystyle{ \pi /2\lt \theta_{1}\lt \theta_{2}\lt {3\pi}/{2}\lt math\gt . Since }[/math]P(\theta,\omega)</math> is negative for [math]\displaystyle{ 0\leq\theta\leq{\pi}/{2} }[/math] or [math]\displaystyle{ {3\pi}/{2}\leq\theta\leq2\pi }[/math] it follows that

[math]\displaystyle{ \frac{\partial P}{\partial\theta}(\theta_{1},0) \gt 0,\quad \mathrm{and}\quad \frac{\partial P}{\partial\theta}(\theta_{2},0) \lt 0.\,\,\,(9) }[/math]

We can determine the type of the equilibrium points by considering the Jacobian matrix, which is given by

[math]\displaystyle{ \left( \begin{matrix} -2\sigma\left\vert \omega\sin\omega-1/\omega\right\vert & \frac{\partial P}{\partial\theta}\\ \omega^{2} & 0 \end{matrix} \right) . }[/math]

Its eigenvalues are

[math]\displaystyle{ \lambda_{1,2}=-\sigma\left\vert \omega\sin\omega-1/\omega\right\vert \pm \sqrt{\sigma^{2}\left\vert \omega\sin\omega-1/\omega\right\vert ^{2} +\frac{\partial P}{\partial\theta}\frac{\omega^{2}}{\sqrt{1+\omega^{4}\cos ^{2}\theta}}}. }[/math]

It follows from (9) that at [math]\displaystyle{ \theta_{1} }[/math] there is one eigenvalue with positive real part and one with negative real part. At [math]\displaystyle{ \theta_{2} }[/math] both eigenvalues have a negative real part. Applying the equilibrium point classification theorems, [math]\displaystyle{ \theta_{1} }[/math] is a saddle point and [math]\displaystyle{ \theta_{2} }[/math] is an attracting node or a spiral.

In summary we have shown that all solutions to equation (7) and hence to equation (6) must tend to either an equilibrium point (at which the velocity of the body is given by the wave phase speed) or to an encircling limit cycle. There can exist at most one encircling limit cycle. The equilibrium points exist only if the drag [math]\displaystyle{ \sigma }[/math] is sufficiently small and they come in an attracting node and saddle pair (a saddle-node bifurcation).