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 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{d\bar{V}_{\tau}}{dt}=-mg\frac{\bar{\eta}^{\prime}}{\sqrt{1+\bar{\eta }^{\prime}{}^{2}}} + F_d\,\,\,(1) (marchenko_original) }[/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{d\bar{V}_{\tau}}{dt} & =\frac{d}{dt}\left( \bar{V}_{\bar{x}} \sqrt{1+\bar{\eta}^{\prime}{}^{2}}\right) \\ & =\frac{d\bar{V}_{\bar{x}}}{dt}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}+\bar {V}_{\bar{x}}\frac{d}{dt}\sqrt{1+\bar{\eta}^{\prime}{}^{2}}\\ & =\frac{d\bar{V}_{\bar{x}}}{dt}\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{d\bar{x}}{dt}. \end{matrix} }[/math]

Substituting this in equation ((marchenko_original)) we obtain

[math]\displaystyle{ m\frac{d\bar{V}_{\bar{x}}}{dt}=-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) (marchenko_travelling) }[/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{dV_{x}}{dt}=-m(g+\left( V_{x}\right) ^{2} \eta^{\prime\prime})\frac{\eta^{\prime}}{1+\eta^{\prime}{}^{2}} + F_d,,\,\,(3) (rumer_correct_added_mass) }[/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{dV}{dt}=-mg\eta^{\prime}+\rho _{w}C_{w}A_{i}\,\,\bigg|V_{w}-V\bigg|\bigg(V_{w}-V\bigg)\,\,\,(4) (linear_force) }[/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) (linear_wave) }[/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)(linear_wave)) into (5)(linear_force)) gives

[math]\displaystyle{ \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) }[/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 ((finalsystemxb)) 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}\,\,\,(6) (autonomous) }[/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}=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 (6)(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

[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 (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

[math]\displaystyle{ \frac{dV}{dt}\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) \ dVd\theta & =\oint_{\Gamma}\left( Pd\theta -QdV\right) \\ & =\int_{0}^{T}\left( P\frac{d\theta}{dt}-Q\frac{dV}{dt}\right) dt\\ & =\int_{0}^{T}\left( PQ-QP\right) dt=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 }\ dVd\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 ((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

[math]\displaystyle{ \begin{matrix} \frac{dV}{dt} & =0,\\ \frac{d\theta}{dt} & =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 (equipoint) }[/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,\\ \frac{\partial P}{\partial\theta}(\theta_{2},0)\ & \lt \ 0. (equiclassi) }[/math]

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

[math]\displaystyle{ \left( \begin{matrix} [c]{cc} -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 ((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

[math]\displaystyle{ \sigma =\frac{\rho_w C_w H}{2\rho_i D}. }[/math]

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,

1. 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:

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

[math]\displaystyle{ \begin{matrix} (zerosof)g(v_{0}) = v_{0} - v(0;2\pi,v_{0}) = v_{0} - v_{1}, \end{matrix} }[/math]

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

[math]\displaystyle{ \begin{matrix} v_{d} = c - \frac{L}{\tilde{t}}. \end{matrix} }[/math]

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}