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