Wave Forcing of Small Bodies
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 begin by comparing the models for wave forcing of small floating bodies which have been derived by Rumer et. al. 1979 and Marchenko 1999. The models 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 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 fashion, the force being due to the difference between the body and water 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]\displaystyle{ x }[/math] direction. In the Marchenko model the coordinate system travels with the wave and the velocity in the Marchenko model is in the tangential direction. A further difference is that Rumer 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]\displaystyle{ 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 system is
where [math]\displaystyle{ \bar{V}_{\tau} }[/math] is the tangential velocity. To compare with equation ((rumer_original)) we need to transform the velocity from the tangential to the [math]\displaystyle{ \bar{x} }[/math] direction. The velocity in the [math]\displaystyle{ \bar{x} }[/math] direction is given by
Therefore
Substituting this in equation ((marchenko_original)) we obtain
which is the Marchenko model with out drag in the [math]\displaystyle{ \bar{x} }[/math] direction. It we 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
It is worth presenting here the corrected Rumer Model in the non-moving coordinate system with the added mass included.
This equation can be applied in situations in which the slope profile is not given by a simple progressing wave.
\bigskip 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 between Rumer and Marchenko about how this force should 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 small body. Marchenko using the velocity in the tangential dircection, while Rumer uses the velocity in the [math]\displaystyle{ x- }[/math]direction. The model for the drag force is therefore given by
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: Theoretical Results
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 theory Shen and Ackley 1991 and Shen and Zhong 2001 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
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
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
where we have assumed that the water depth in infinite. Substituting equation ((linear_wave)) into ((linear_force)) gives \begin{gather} \left( 1+C_{m}\right) m\frac{dV}{dt}\,=-m\,g\,\frac{kA}{2}\cos(kx-\omega t) (finalsystemxb)\\ + \,\rho_{w}AC_{w}|kcA\sin(k(x-ct))-V|(kcA\sin(k(x-ct))-V)\nonumber \end{gather}
Non-dimensionalisation
We will now non-dimensionalised 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 ((finalsystemxb)) becomes the following system
of equations
where the variables are now non-dimensional and where
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
Equation ((autonomous)) 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
It should also be noted that all the constants must be positive.
The behaviour of the solutions
Analytic solutions for system (linear_force) 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{ dV/dt }[/math] can therefore be approximated by
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, perko91, 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.
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
which is known as Bendixon's criterion. However,
which is always negative except on the line [math]\displaystyle{ V=\omega\sin{\theta} }[/math]. Therefore
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 ((autonomous)) and hence to equation ((finalsystemxb)) 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
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.
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
(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
We can determine the type of the equilibrium points by considering the Jacobian matrix, which is given by
Its eigenvalues are
It follows from ((equiclassi)) 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 (perko91), [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 ((autonomous)) and hence to equation ((finalsystemxb)) 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).
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]\displaystyle{ \omega }[/math] [math]\displaystyle{ \tau }[/math], and [math]\displaystyle{ \sigma\lt math\gt , which have to be determined. }[/math]\omega</math> can be calculated directly from the wavelength and the amplitude of the wave by [math]\displaystyle{ \omega= \sqrt{2 \pi A/L } }[/math]. To determine [math]\displaystyle{ \sigma }[/math], the relative density of the ice to that of the water, the non-dimensionalised height [math]\displaystyle{ 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]\displaystyle{ C_m=0.08 }[/math] so [math]\displaystyle{ \tau=1/(1.08)=0.93 }[/math]. They assumed that the floating body was of constant cross section so that the mass was [math]\displaystyle{ m=\rho_i A D }[/math] where [math]\displaystyle{ \rho_i }[/math] is the density and [math]\displaystyle{ D }[/math] is the body thickness. This means that we can express [math]\displaystyle{ \sigma }[/math] as
Shen \& Zhong further assumed that [math]\displaystyle{ H/\lambda = 0.02 }[/math], that [math]\displaystyle{ C_w\lambda/D = 0.1 }[/math] where [math]\displaystyle{ \lambda }[/math] is the wavelength and that [math]\displaystyle{ \rho_w/\rho_i = 1/0.9 }[/math] This gives us a value for [math]\displaystyle{ \sigma }[/math] of [math]\displaystyle{ \sigma = 0.02/1.8=0.011 }[/math]. The value of [math]\displaystyle{ \omega }[/math] is found from [math]\displaystyle{ \omega = \sqrt{\pi H/\lambda} }[/math] so that [math]\displaystyle{ \omega = \sqrt{2\pi\times 0.01} = 0.25 }[/math]. We can see that [math]\displaystyle{ \omega^2/(\omega-1/\omega)^2 = 0.0045 }[/math] which is less than [math]\displaystyle{ \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]\displaystyle{ \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]\displaystyle{ \omega }[/math] is [math]\displaystyle{ 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]\displaystyle{ \omega }[/math] is given by 0.8. In this case, two equilibrium points, [math]\displaystyle{ (\theta_{1}, {1}/{\omega}) }[/math] and [math]\displaystyle{ (\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]\displaystyle{ -1 }[/math]. [math]\displaystyle{ \theta_{1} }[/math] and [math]\displaystyle{ \theta_{2} }[/math] can be found by solving equation ((equipoint)), for example by finding the zeros of [math]\displaystyle{ P(\theta,{1}/{\omega }) }[/math] numerically with Newton's method. The starting points can be chosen as [math]\displaystyle{ x_{0}={\pi}/{2} }[/math] and as [math]\displaystyle{ 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]\displaystyle{ V(\theta;\theta_{0},v_{0}) }[/math] denotes the velocity [math]\displaystyle{ V(\theta) }[/math] satisfying the initial condition [math]\displaystyle{ v(\theta_{0})=v_{0} }[/math]. It should be noted that the dependent variable [math]\displaystyle{ 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]\displaystyle{ v_{0} }[/math] at [math]\displaystyle{ \theta_{0} = 2\pi }[/math] such that
the velocity at [math]\displaystyle{ \theta= 0 }[/math], [math]\displaystyle{ v(0;2\pi,v_{0}) = v_{1} }[/math], is equal to the initial velocity, [math]\displaystyle{ v_{0} = v_{1} }[/math]? This involves the multiple solving of the sub-problem:
- Given an initial condition [math]\displaystyle{ (2\pi, v_{0}) }[/math], what is the velocity [math]\displaystyle{ 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]\displaystyle{ \theta=0 }[/math]: The closest negative value, [math]\displaystyle{ \theta_{-} }[/math], and the closest positive
one, [math]\displaystyle{ \theta_{+} }[/math]. Using these values, a good approximation of the velocity at
[math]\displaystyle{ \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]\displaystyle{ (2\pi, v_{0}) }[/math]. The velocity [math]\displaystyle{ 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]\displaystyle{ v_{1} }[/math] is used as the new initial condition: The system is solved again with initial condition [math]\displaystyle{ (2\pi, v_{1}) }[/math] and the whole method is iterated. The iteration stops when [math]\displaystyle{ |v_{n} - v_{n-1}| }[/math] is smaller than a tolerance [math]\displaystyle{ \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]\displaystyle{ v_{0} = v_{1} }[/math] can be found by finding the zero of
where [math]\displaystyle{ v_{0} }[/math] is the initial velocity at [math]\displaystyle{ \theta_{0}=2\pi }[/math] and [math]\displaystyle{ v_{1} }[/math] is the corresponding velocity at [math]\displaystyle{ \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]\displaystyle{ \tilde{t} }[/math] the ice floe needs to travel from one wavecrest to the next wavecrest using the formula
However, a small change in the value of [math]\displaystyle{ \tilde{t} }[/math] in this calculation results in a great variation in [math]\displaystyle{ 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]\displaystyle{ \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).
