Difference between revisions of "Wave Forcing of Small Bodies"

From WikiWaves
Jump to navigationJump to search
 
(22 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 begin by comparing the models for wave forcing of small floating bodies
+
We present the models for wave forcing of small floating bodies
which have been derived by [[Rumer et. al. 1979]] and [[Marchenko 1999]]. The models
+
which have been derived by [[Marchenko 1999]]. Marchenko decomposed the force acting of the small
are similar, both Rumer and Marchenko decomposed the force acting of the small
 
 
body into two components due to the drag force between the water and the body
 
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.
 
and the gravity force due to the body sliding down the surface of the wave.
The drag force was modelled by Rumer and Marchenko in almost indentical
+
The drag force is due to the difference between the body and water
fashion, the force being due to the difference between the body and water
+
velocities squared. In the Marchenko model the
velocities squared. However in the gravity sliding force Rumer and Marchenko
 
differ in their derivation. In the Rumer model the coordinate system is fixed
 
and the velocity is given in the <math>x</math> direction. In the Marchenko model the
 
 
coordinate system travels with the wave and the velocity in the Marchenko
 
coordinate system travels with the wave and the velocity in the Marchenko
model is in the tangential direction. A further difference is that Rumer
+
model is in the tangential direction.  
includes a term for the added mass. For the purpose of comparision, we will
 
convert both equations to a coordinate system which is travelling with the
 
wave, velocity will be given in the <math>x</math> direction and the added mass will be
 
neglected. We will also consider the models without the drag force so that we
 
concentrate on the sliding force. 
 
 
 
 
 
 
 
==The Marchenko Model without Drag==
 
  
 
The model for the sliding force given by Marchenko in the moving co-ordinate
 
The model for the sliding force given by Marchenko in the moving co-ordinate
 
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}}} (marchenko_original)
+
}^{\prime}{}^{2}}} + F_d\,\,\,(1)  
 
</math></center>
 
</math></center>
where <math>\bar{V}_{\tau}</math> is the tangential velocity. To compare with equation
+
where <math>\bar{V}_{\tau}</math> is the tangential velocity and
((rumer_original)) we need to transform the velocity from the tangential
+
<math>F_d</math> is the drag force (which we will introduce later).  
to the <math>\bar{x}</math> direction. The velocity in the <math>\bar{x}</math> direction is given
+
The velocity in the <math>\bar{x}</math> direction is given
 
by
 
by
 
<center><math>
 
<center><math>
Line 78: 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 ((marchenko_original)) 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}}. (marchenko_travelling)
+
{}^{2}} + F_d.\,\,\,(2)  
 
</math></center>
 
</math></center>
which is the Marchenko model with out drag in the <math>\bar{x}</math> direction. It we
+
which is the Marchenko model in the <math>\bar{x}</math> direction.  
comparing equations ((rumer_travelling)) and equation
 
((marchenko_travelling)) we see that Markenko and Rumers models are not
 
the same. In the folowing section we will derive the equations of motion using
 
Hamilton's principle to determine which of the equations is correct.
 
 
 
 
 
 
 
  
  
==Corrected Rumer Model in the non-moving frame with added mass==
+
===Including added mass===
  
It is worth presenting here the corrected Rumer Model in the non-moving
+
[[Rumer et. al. 1979]] included added mass and this makes the model
coordinate system with the added mass included.
 
 
<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}}
+
\eta^{\prime\prime})\frac{\eta^{\prime}}{1+\eta^{\prime}{}^{2}} + F_d,\,\,\,(3)
(rumer_correct_added_mass)
 
 
</math></center>
 
</math></center>
This equation can be applied in situations in which the slope profile is not
 
given by a simple progressing wave.
 
  
==\bigskip 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
 
force is the more difficult force to model because it depends on an unknown
 
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
 
factor which models the friction force between the body and the water. However
there is a consensus between Rumer and Marchenko about how this force should
+
there is a consensus that the drag force should be proprotional to the square
be modelled.  Both model the drag force as a being proprotional to the square
 
 
of the velocity difference of the water particles at the water surface and the
 
of the velocity difference of the water particles at the water surface and the
small body. Marchenko using the velocity in the tangential dircection, while
+
small body. The model for the drag force is
Rumer uses the velocity in the <math>x-</math>direction. The model for the drag force is
 
 
therefore given by
 
therefore given by
 
<center><math>
 
<center><math>
Line 134: Line 108:
 
co-ordiante system.
 
co-ordiante system.
  
=The Drift Velocity: Theoretical Results=
+
==The Drift Velocity for a Linear Sinusoidal Wave==
 
 
One of the priciple aims of deriving the system of equations for the wave
 
forcing of a small floating body was the determination of the drift velocity.
 
Calculating the drift velocity was the principle aim of the papers
 
[[Shen and Ackley 1991]] and [[Shen and Zhong 2001]]. However these papers claimed that
 
the drift velocity was a function of the initial conditions for the case of a
 
regular sinusoidal wave which we will show here is false.
 
 
 
==The systems of equations 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
theory [[Shen and Ackley 1991]] and [[Shen and Zhong 2001]] then it makes sense to work
+
theory then it makes sense to work
 
derive the system of equations under the same assumptions which underlie the
 
derive the system of equations under the same assumptions which underlie the
 
linear wave theory. This means that the tangential and <math>x</math> directions should
 
linear wave theory. This means that the tangential and <math>x</math> directions should
Line 153: 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) (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 160: 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)), (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 168: 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
((linear_wave)) into ((linear_force)) gives  
+
(4) into (5) gives  
\begin{gather}
+
<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)\nonumber
+
+ \,\rho_{w}AC_{w}|kcA\sin(k(x-ct))-V|(kcA\sin(k(x-ct))-V). \,\,\,(6)
\end{gather}
+
</math></center>
  
  
==Non-dimensionalisation==
+
===Non-dimensionalisation===
  
We will now non-dimensionalised equation ((finalsystemxb)). This will
+
We will now non-dimensionalise equation ((finalsystemxb)). This will
 
serve two purposes, the first to simplify the equations and the second to
 
serve two purposes, the first to simplify the equations and the second to
 
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
+
<math>\sqrt{{H}/{2g}}</math> so that equation (6) becomes the following system
<center><math>
 
\tilde{x} = \frac{2x}{H}, \mbox{\hspace{1cm}}
 
\tilde{t} = \sqrt{\frac{2g}{H}}t, \mbox{\hspace{1cm}}
 
\rm{and} \mbox{\hspace{1cm}}
 
\tilde{V} = \sqrt{\frac{1}{2gH}}V \,
 
</math></center>
 
Equation ((finalsystemxb)) becomes the following system
 
 
of equations
 
of equations
 
<center><math>
 
<center><math>
 
\begin{matrix}
 
\begin{matrix}
[c]{c}
 
 
\frac{d\tilde{V}}{d\tilde{t}}\,=-\tau\omega^{2}\,\cos{\theta}+\sigma\tau\left\vert \omega
 
\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)  \\
 
\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}
+
\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 213: Line 169:
 
<center><math>
 
<center><math>
 
\begin{matrix}
 
\begin{matrix}
[c]{c}
 
 
P(V,\theta)=-\omega^{2}\,\cos{\theta}+\sigma\left\vert \omega\sin{\theta
 
P(V,\theta)=-\omega^{2}\,\cos{\theta}+\sigma\left\vert \omega\sin{\theta
 
}-V\right\vert \,\left(  \omega\sin{\theta}-V\right)  \\
 
}-V\right\vert \,\left(  \omega\sin{\theta}-V\right)  \\
Line 219: Line 174:
 
\end{matrix}
 
\end{matrix}
 
</math></center>
 
</math></center>
Equation ((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 230: 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 (linear_force) cannot be obtained easily.
+
Analytic solutions for system (4) cannot be obtained easily.
 
Starting with the simplified [[Rumer et. al. 1979]] model (without centripetal force),
 
Starting with the simplified [[Rumer et. al. 1979]] model (without centripetal force),
 
[[Shen and Zhong 2001]] found approximate analytic solutions for a few special
 
[[Shen and Zhong 2001]] found approximate analytic solutions for a few special
Line 238: 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 247: Line 202:
 
lives on the cylinder -<math>\pi\leq\theta\leq\pi</math>, <math>-\infty<V<\infty,</math> cannot
 
lives on the cylinder -<math>\pi\leq\theta\leq\pi</math>, <math>-\infty<V<\infty,</math> cannot
 
leave a bounded region of the cylinder.  This means that we can use the
 
leave a bounded region of the cylinder.  This means that we can use the
Poincare-Bendixson theorem, [[perko91]], which in the case of a cylinder
+
[[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
 
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
 
equilibrium point. It should be noted that there are two possible limit cycles
Line 259: Line 214:
 
by <math>\Gamma,</math> i.e.
 
by <math>\Gamma,</math> i.e.
 
<center><math>
 
<center><math>
\Gamma:=\left\{  \binom{V(t)}{\theta(t)}\,\big|\,0\leq t\leq T\right\}
+
\Gamma=\left\{  ({V(t)},{\theta(t)})\,\big|\,0\leq t\leq T\right\}
 
</math></center>
 
</math></center>
 
and <math>\Gamma\left(  0\right)  =\Gamma\left(  T\right)  </math> in a simply connected
 
and <math>\Gamma\left(  0\right)  =\Gamma\left(  T\right)  </math> in a simply connected
Line 265: 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 278: 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 290: 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 298: 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 306: 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 330: 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 341: 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 355: 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
 
cycle. The equilibrium points exist only if the drag <math>\sigma</math> is sufficiently
 
cycle. The equilibrium points exist only if the drag <math>\sigma</math> is sufficiently
 
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}
 
 
 
 
 
=Summary and discussion=
 
 
 
In this work, the motion of a small body under the influence of water waves
 
was investigated. Slope sliding models and the linear water wave theory were
 
used to describe the motion of a small ice floe floating on a periodic water wave.
 
 
 
To model how the body is affected by the wave field, two similar slope sliding
 
models were considered: The system of equations derived by [[Rumer et. al. 1979]] uses
 
velocities in horizontal direction and a stationary coordinate system while
 
the system of equations for the motion given by [[Marchenko 1999]] is
 
formulated for the velocities in the tangential direction measured in a
 
coordinate system moving with the same speed as the wave. Both systems were
 
derived explicitly. In a thorough comparison, also including the derivation of
 
a system of equations using Hamilton's principle, it was shown that the two
 
systems are not equivalent since the [[Rumer et. al. 1979]] model does not include the
 
centripetal force correctly.
 
 
 
This error occurred since in the system of equations a projection of the
 
inertia in vertical direction onto the horizontal direction was used, but the
 
centripetal force due to the curvature of the wave was ignored. To correct the
 
system of equations, there are two options: Either all forces introduced by
 
[[Rumer et. al. 1979]], including the inertia in vertical and horizontal direction,
 
are projected onto the tangential instead of the normal (in which case the
 
centripetal force has no influence since it points perpendicular to the
 
tangential), or the centripetal force is included and all forces are projected
 
onto the horizontal direction (in which case the vertical inertia has no
 
influence). Both options lead to the same system of equations as the one
 
derived by using Hamilton's principle or [[Marchenko 1999]]'s model.
 
 
 
[[Marchenko 1999]]'s system of equations and its equivalent formulated for the
 
horizontal velocity in a stationary reference frame were applied to a small
 
ice floe floating on periodic water waves. The results from the linearised
 
water theory summarised in the first chapter were substituted in the system of
 
equations. It was shown that under the assumption of a long wave with very
 
small amplitude the centripetal term as well as parts of the drag force can be
 
neglected, which results in a system which is the same as the simplified
 
[[Rumer et. al. 1979]] model. Therefore this simplified [[Rumer et. al. 1979]] model, which
 
has been used by [[Shen and Ackley 1991]], [[frankshen93]], [[meylanphd94]]
 
and others to model the movement of ice floes, is a good approximation as long
 
as the wave slope it not too steep.
 
 
 
As suggested by [[meylanphd94]], the general behaviour of solution was
 
investigated. It was shown that after sufficiently long time the ice floe
 
either travels with the same velocity as the wave, is therefore captured by
 
the wave, or is moving from one wave crest to the next in always the same way.
 
For a few settings numerical solutions were presented. The solutions obtained
 
from the improved system of equations were compared with solutions obtained
 
from the simplified [[Rumer et. al. 1979]] model without centripetal force.
 
 
 
A numerical procedure to calculate the limit cycle was described. This
 
procedure was used to determine the long term drift velocity. It was shown
 
that increasing the drag coefficient leads to a higher drift velocity. The
 
same result can be obtained by decreasing the height of the ice floe. Without
 
the drag force, the long term drift velocity is not unique anymore and only
 
depends on the initial ice floe velocity.
 
 
 
The solutions obtained from the [[Rumer et. al. 1979]] model are very similar to those
 
of the corrected model if the amplitude of the wave is very small compared to
 
the wavelength, as required for the linear water wave theory. For larger
 
amplitudes the solutions begin to differ since then the influence of the
 
centripetal force and the difference in the drag forces increases (due to the
 
steeper slope and higher curvature of the wave surface).
 
 
 
[[Category:Linear Water-Wave Theory]]
 

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).