Ulrich Razafison
email: ulrich.razafison@math.cnrs.fr

Stephane Cordier
email: stephane.cordier@univ-orleans.fr

Olivier Delestre
email: olivierdelestre41@yahoo.fr

Frédéric Darboux
email: frederic.darboux@orleans.inra.fr

Carine Lucas
email: carine.lucas@univ-orleans.fr

Francois James
email: francois.james@univ-orleans.fr

Stéphane Cordier

François James

A shallow water

Keywords: Overland flow, Shallow water equations, Furrows, Friction We introduce AMS Classification: 93A30, 81T80, 58J45, 35L65

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

During rainfalls, overland flow on cultivated lands induces problems at watershed scale for soil conservation (decrease of soil thickness by erosion, nutrient losses), infrastructures (flooding and destruction of roads and buildings), preservation of water quality (drinking water) and sustainability of aquatic ecosystems (chemical pollution).

These troubles can be prevented by improving watershed management in connection with overland flow. This requires to simulate well the water flux at the outlet but also to have a good prediction of the spatial distribution of the water flux and velocity over the whole watershed. However, the current hydrological models have a low efficiency in predicting overland flow within small watersheds (see [START_REF] Beven | On undermining the science?[END_REF][START_REF] Jetten | Evaluation of field-scale and catchmentscale soil erosion models[END_REF][START_REF] Jetten | Erosion models: quality of spatial predictions[END_REF]). In agricultural watersheds, one of the main difficulties is that flow directions are controlled not only by the topography, but also, by ditches along the field boundaries, and, inside the fields, by ridges and furrows created by tillage operations. The flow pattern is clearly the result of the interaction between these objects [START_REF] Souchère | Effects of tillage on runoff directions: consequences on runoff contributing areas within agricultural catchments[END_REF], but the way they interact remains mostly unspecified. Therefore, one must improve the understanding of this interaction in order to better predict the spatial and temporal distribution of overland flow and so to improve the decisions made by watershed managers.

In this paper we focus on the interaction between topography and furrows, a feature encountered in almost all cultivated lands. This interaction can be seen as the interaction between three types of roughness. The topography is the roughness of the Earth and is described in Digital Elevation Maps with a horizontal resolution larger than one meter, and commonly of ten meters and more. The furrows are the roughness due to agricultural practices and create a strong directional heterogeneity inside a field. They are characterized by their wavelength (of about one to a few decimeters), their amplitude (of a few centimeters to one decimeter) and their direction. Finally the random roughness, due to soil aggregates and clods, is homogeneous in space and has an amplitude of a few millimeters to about one decimeter. To our knowledge, most of the works on the interaction between roughness and flow have been dedicated to topography (see [START_REF] Quinn | The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models[END_REF][START_REF] Tarboton | A new method for the determination of flow directions and upslope areas in grid digital elevation models[END_REF]) or to random roughness (see [START_REF] Darboux | Effect of depression storage capacity on overland-flow generation for rough horizontal surfaces: Water transfer distance and scalling[END_REF][START_REF] Jin | Estimating Manning's roughness coefficient for shallow overland flow in non-submerged vegetative filter strips[END_REF][START_REF] Gilley | Darcy-Weisbach roughness coefficients for overland flow[END_REF]).

Few works are dealing with furrows, and among them, most are concerned with the storage capacity of the furrows, i.e. the amount of water stored in the puddles created by the furrows (for instance [START_REF] Planchon | Microrelief induced by tillage: measurement and modelling of Surface Storage Capacity[END_REF]). These works do not consider the water flowing on the soil surfaces but only the water stored in puddles. The few works considering both overland flow and the furrows-topography interaction are empirical studies [START_REF] Souchère | Effects of tillage on runoff directions: consequences on runoff contributing areas within agricultural catchments[END_REF][START_REF] Takken | Effect of tillage and runoff and erosion patterns[END_REF]. They lead to empirical laws giving an on/off prediction: the predicted flow direction is either the direction of the topographic slope or the furrow direction, while in reality water can flow in both directions at the same time. Moreover, these laws are limited by their empirical basis.

To be of practical use, a model accounting for the effects of furrows on overland flow direction must not require an explicit representation of the furrows: that would imply the use of a digital topographic map with a horizontal resolution of about a centimeter for the whole watershed, which already covers an area of about one square kilometer for a small one. Such digital maps are not available and, even if available, will require too much computation resources. of the furrows on overland flow. Numerical results are presented. The model is a first step in an attempt to predict overland flow directions controlled by furrows and topography without representing the furrows explicitly. Indeed, the furrows are known only trough their average amplitude, wavelength and direction. In this paper, the furrow direction is kept perpendicular to the slope. Our model is based on the Shallow Water equations that are widely used to describe flows in rivers, in ocean and overland flow among other applications.

The outline of the paper is as follows. In the next section, we present first the Shallow Water model. Then we propose a new model where we add a new friction term to account for the effects of the furrows on overland flow. Section 3 describes the numerical scheme used to solve the model, and, in section 4, we present and discuss the numerical results that we obtain with our model. Conclusions are outlined in Section 5.

The mathematical models

The starting point is the 2D classic Shallow Water system [START_REF] De Saint-Venant | Théorie du mouvement non-permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF] in a bounded domain Ω: (2.1)

                 ∂h ∂t + ∂(hu) ∂x + ∂(hv) ∂y = R, ∂ (
For t > 0 and x = (x, y) ∈ Ω, the unknowns are the water height h = h(t, x) and the horizontal flow velocity u = u(t, x) = (u(t, x), v(t, x)) T . Furthermore, Z(x) describes the bottom topography of the domain and therefore h + Z is the level of the water surface (see Figure 1). In equations (2.1), g is the acceleration due to the gravity and R is the rainfall intensity. We refer for instance to [START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF][START_REF] Ferrari | A new two-dimensional shallow water model including pressure effects and slow varying bottom topography[END_REF][START_REF] Marche | Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects[END_REF] for a derivation of the Shallow Water system originating from the free surface Navier-Stokes equations.

For the friction term, we choose the Manning law with k the Manning coefficient. We also denote q(t, x) = (q x (t, x), q y (t, x)) T = h(t, x)u(t, x) the water flux. Now we consider a rectangular domain Ω = ℓ × L and a topography Z with furrows.

We suppose that the topography is an inclined plane with sinus furrows, and that the geometry of the furrows is known through their amplitude and their wavelength. Note that realistic (measured) furrows will be a slightly different shape due the existence of random roughness. However, random roughness, being isotropic, does not affect flow direction at the scale of the furrows. We also suppose that the furrows are perpendicular 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 1. The direction of the flow is parallel to the length of the domain Ω with respect to y (pseudo-1D case) and consequently perpendicular to the furrows.

11111111111111111111111111111111111111
2. We only consider fluvial flows which means that |u| < √ gh.

3. Infiltration and soil erosion are not taken into account.

Under such assumptions, the furrows overflow at the same time during rainfall events or one after the other in the case of an inflow from upstream.

We aim at proposing a model that takes into account the effects of the furrows without explicitly representing them in the topography Z. In other words, we want to find an equivalent model to the Shallow Water system on Ω, that would be used at a macroscopic scale, i.e. on a topography which is only an inclined plane. We want to force the flow to slow down when its depth is smaller than the value corresponding to the water height that can be trapped in the furrows. The idea of this article is to model this effect caused by the furrows through an additional friction term that forces the flow to slow down for small water depth. To that end, we first introduce h F the average height of the water trapped in the furrows. This value is given by

h F = V /(L F × ℓ) [m] (2.2)
where V is the volume of trapped water in a furrow, L F its wavelength (see Figure 3) and ℓ is the length of the domain Ω (with respect to x). Note that the value of h F only depends on three parameters: the slope of the domain, the furrows average amplitude and the furrows average wavelength. 00000000000000000 00000000000000000 00000000000000000 11111111111111111 11111111111111111 11111111111111111

y 0 + L F z y 0 Z (y) h(t, y) y
Trapped water (volume V ) Next we consider the following additional friction coefficient:

K(h) = K 0 exp -h + h F C h F , (2.3)
where C is a characteristic constant, increasing function of the small random variations of the height of the furrows and K 0 is a coefficient we determine in the following.

In Figure 4, the general shape of K(h) is plotted for h F = 0.01 m. We clearly see that K(h) is large for h ≤ h F . This shows that when the water height h is lower than the average height of the furrows h F then, thanks to K(h), the flow is slowed down.

Remark 2.1. Let us explain in a few words relation (2.3).

1. If h F tends to 0, then K(h) also tends to 0 for any h > 0. In other words, the additional friction coefficient disappears when there are no furrows. If C tends to 0, then we obtain the empirical models that are usually used. These models consist in giving an on/off prediction of the furrows-topography interaction (see [START_REF] Souchère | Effects of tillage on runoff directions: consequences on runoff contributing areas within agricultural catchments[END_REF][START_REF] Takken | Effect of tillage and runoff and erosion patterns[END_REF]); more precisely, while the critical water height is not attained, there is no flow, and after this threshold, the furrows are not taken into account.

The new Shallow Water model we introduce here (see (2.4)) can be seen as an improvement of these models.

Finally, we propose the following new Shallow Water model with a "furrows-friction"

coefficient:                  ∂h ∂t + ∂(hu) ∂x + ∂(hv) ∂y = R, ∂(hu) ∂t + ∂(hu 2 ) ∂x + ∂(huv) ∂y + gh ∂h ∂x + gh ∂Z ∂x + gk 2 h -1/3 |u|u = 0, ∂(hv) ∂t + ∂(huv) ∂x + ∂(hv 2 ) ∂y + gh ∂h ∂y + gh ∂Z ∂y + gk 2 h -1/3 |u|v + K(h)hv = 0. (2.4) Remark 2.2.
1. Note that, since the furrows are perpendicular to the slope, the additional friction law K(h)hv only appears in the third equation of (2.4) and therefore, it only acts on the flow in the y-axis direction. This assumption is not restrictive: in general, the direction of the furrows is constant on each agricultural field, and, if necessary, we apply a rotation to get the equations for an arbitrary direction of the furrows.

2. The form of the new friction law is arbitrarily chosen. The general form of friction laws is Kh α |u| β u where α and β are positive real numbers. For example we get Manning's law for (α, β) = (-1/3, 1) and Darcy-Weisbach's law for (α, β) = (1, 1).

Note that these laws are empirical, are obtained considering stationary flows and that their validity is still discussed among hydrologists (e.g. [START_REF] Smith | Applying flow resistance equations to overland flows[END_REF]).

For the numerical experiments presented in section 4 with system (2.4), we chose (α, β) = (1, 0) but, we could have made another choice, should we change the value of K 0 or the form of K(h).

At this point, let us mention that, since shallow water flows can also be described by the so-called multi-layer Shallow Water system (see for instance to [START_REF] Audusse | A multilayer Saint-Venant model: derivation and numerical validation[END_REF][START_REF] Audusse | Finite-volume solvers for a multilayer Saint-Venant system[END_REF][START_REF] Castro | A Q-scheme for a class of systems of coupled conservation laws with source term. Application to a two-layer 1-D shallow water system[END_REF]] for a derivation and numerical studies), we can propose another approach based on multi-layer models to take into account the effects of the furrows on overland flows.

In this work, we introduce the following two-layer like model:

     if h(t, x) ≤ h F , then u(t, x) = 0 and h(t, x) = Rt, if h(t, x) > h F , then solve (2.1
) with an inclined plane topography.

(

In (2.5), the lower layer corresponds to the filling up of the furrows; note that the upper layer is active only when the furrows overflow. The initial conditions for the upper layer are then u(0, x) = 0 and h(0, x) = ĥ -h F , where ĥ is the water height at the overflow time. Note that, in one dimension, this model provides more satisfactory results than the model (2.4) (see section 4). But its extension to more complex two-dimensional problems requires a careful modelling of the coupling between the two layers and it is more difficult than the extension of the model (2.4).

Numerical results

In this section, we explain the numerical scheme we used in our numerical simulations.

The Shallow Water system is not so easy to solve. In hydrology, Mac Cormack scheme is usually used for overland flow simulation (see among others [START_REF] Esteves | Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values[END_REF][START_REF] Fiedler | A numerical method for simulating discontinuous shallow flow over an infiltrating surface[END_REF]). But it is not well adapted to this system because of several problems, such as, for example, the preservation of the positivity of the water height and of the steady states or the behavior at the wet/dry interface. To that end, we used well-balanced schemes (so-called since [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equation[END_REF]) based on the hydrostatic reconstruction [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. This finite volume scheme has shown to be adapted to overland flow simulation at small scales [START_REF] Delestre | Simulation of rain-water overlandflow[END_REF][START_REF] Delestre | Simulation of rainfall events and overland flow[END_REF][START_REF] Delestre | Rain water overland flow on agricultural fields simulation[END_REF].

To make this presentation simpler, we describe the numerical scheme on the classical the one-dimensional Shallow Water model with variable topography and Manning's friction law :

       ∂h ∂t + ∂(hu) ∂x = R ∂(hu) ∂t + ∂ ∂x hu 2 + g h 2 2 = -gh ∂Z ∂x -gk 2 h -1/3 |u|u. (3.6)
The model (3.6) can be written into a conservative form

∂U ∂t + ∂F (U ) ∂x = S 0 (U ) + S f (U ), (3.7)
where

U = h hu = h q , F (U ) = hu hu 2 + g h 2 2 , S 0 (U ) = R -gh ∂Z ∂x and S f (U ) = 0 -gk 2 h -1/3 |u|u .
System (3.7) is discretized using finite volume method for hyperbolic conservation laws.

We introduce a space-time grid where the space and the time steps are respectively ∆x and ∆t. We set x i = i∆x, t n = n∆t and C i = x i-1/2 , x i+1/2 . We denote by U n i the approximation of the average of U (t n , x) over the cell C i , namely,

U n i ≃ 1 ∆x C i U (t n , x)dx.
Considering for the moment only the homogeneous part of (3.7), then the finite volume scheme is of the form

U n+1 i -U n i + ∆t ∆x (F n i+1/2 -F n i-1/2 ) = 0,
where

F n i+1/2 = F(U n i , U n i+1 )
is the HLL numerical flux (see for instance [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]) through the interface between C i and C i+1 . Note that the HLL flux is defined by

F(U l , U r ) =              F (U l ) if 0 < c 1 , c 1 F (U l ) -c 2 F (U r ) c 2 -c 1 + c 1 c 2 c 2 -c 1 (U r -U l ) if c 1 < 0 < c 2 , F (U r ) if c 2 < 0,
where c 1 < c 2 are given by c 1 = inf

U =U l ,Ur inf j=1,2 λ j (U ) , c 2 = sup U =U l ,Ur sup j=1,2 λ j (U ) and
where

λ 1 (U ) = u - √ gh, λ 2 (U ) = u + √
gh are the eigenvalues of the Jabobian matrix

F ′ (U ).
In order to have a second order accuracy scheme, we use the modified ENO reconstruction (see [START_REF] Bouchut | Efficient numerical finite volume schemes for shallow water models[END_REF]) defined as follows

h i-1/2+ = h i - ∆x 2 D enom h i , h i+1/2-= h i + ∆x 2 D enom h i , u i-1/2+ = u i - h i+1/2- h i ∆x 2 D enom u i , u i+1/2-= u i + h i-1/2+ h i ∆x 2 D enom u i
with, for a spatially discretized function s,

D enom s i = minmod(D eno s i , 2θ enom D mm s i ) where minmod(x, y) =        min(x, y) if x, y 0 max(x, y) if x, y 0 0 otherwise, D eno s i = minmod s i -s i-1 ∆x + θ eno ∆x 2 D 2 s i-1/2 , s i+1 -s i ∆x -θ eno ∆x 2 D 2 s i+1/2 , D 2 s i+1/2 = minmod s i+1 -2s i + s i-1 ∆x 2 , s i+2 -2s i+1 + s i ∆x 2 , D mm s i = minmod s i -s i-1 ∆x , s i+1 -s i ∆x with θ eno , θ enom ∈ [0, 1]
. Note that for θ eno = 0, this reconstruction is exactly the usual MUSCL reconstruction.

In order to take into account the topography while preserving the steady state of a lake at rest, that is

h + Z = cst and u = 0,
we use the hydrostatic reconstruction developed in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. First we need to define the reconstructed values z i+1/2-and z i-1/2+ that can be deduced from the reconstructed values h i-1/2+ , h i-1/2+ and the following reconstruction of z + h:

(z + h) i-1/2+ = z i + h i - ∆x 2 D enom (z i + h i )
and

(z + h) i+1/2-= z i + h i + ∆x 2 D enom (z i + h i ).
Next the hydrostatic reconstruction consists in defining the following new values:

h i+1/2 l = max(0, h i+1/2-+ z i+1/2--max(z i+1/2-, z i+1/2+ )), h i+1/2 r = max(0, h i+1/2+ + z i+1/2+ -max(z i+1/2-, z i+1/2+ )), U i+1/2 l = h i+1/2 l h i+1/2 l u i+1/2- , U i+1/2 r = h i+1/2 r h i+1/2 r u i+1/2+ .
The positive parts in the definitions of h i+1/2 l and h i+1/2 r (we have max(0, •)) insure the positivity of the water height. Therefore the scheme can be written into the form

U n+1 i -U n i + ∆t ∆x (F n i+1/2 l -F n i-1/2 r -F c n i ) = 0, F n i-1/2 r = F(U n i-1/2 l , U n i-1/2 r ) +   0 g 2 (h n i-1/2+ ) 2 -(h n i-1/2 r ) 2   , F n i+1/2 l = F(U n i+1/2 l , U n i+1/2 r ) +   0 g 2 (h n i+1/2-) 2 -(h n i+1/2 l ) 2   , and
F c n i =   0 - g 2 h n i-1/2+ + h n i+1/2- z n i+1/2--z n i-1/2+   .
The term F c n i is added to obtain a well-balanced and consistent scheme (see [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]). Now in order to have a second order scheme in time, we use the Heun method,

U i n+1 = U n i + ∆t Φ(U n i ) U n+2 i = U n+1 i + ∆t Φ( U n+1 i ) U n+1 i = U n i + U n+2 i 2 , (3.8)
where

Φ(U n i ) = 1 ∆x (F n i+1/2 l -F n i-1/2 r -F c n i ).
Concerning the Manning friction term, we follow [START_REF] Fiedler | A numerical method for simulating discontinuous shallow flow over an infiltrating surface[END_REF][START_REF] Bristeau | Boundary conditions for the shallow water equations solved by kinetics schemes[END_REF] and introduce a semi-implicit treatment of this term. Then the scheme is modified as follows

• Solve the Shallow Water system

U * i = U n i + ∆t Φ(U n i ). • Compute U n+1 i =   h n+1 i h n+1 i u n+1 i 
 by solving the Manning friction term

   h n+1 i h * i u n+1 i -u * i ∆t    = Sf ∆t (U * i ) ≡    h * i -gk 2 q n+1 i |q * i | h n i (h n+1 i ) 4/3    .
• Solve the Shallow Water system

U * * i = U n+1 i + ∆t Φ( U n+1 i ). • Compute U n+2 i =   h n+2 i h n+1 i u n+2 i   by solving the Manning friction term    h n+2 i h * * i u n+2 i -u * * i ∆t    = Sf ∆t (U * * i ).
consider two types of test cases: in the first one, we only take into account rainfalls and in the second one the water comes from upstream. For these two numerical experiments, the "exact" (or reference) solution is the one of the Shallow Water system (2.1) with a precise description of a topography with furrows. The domain Ω we consider here is Ω = ℓ × L, where ℓ = 0.2 m and L = 4 m (see Figure 2). We assume that the plane topography has a constant slope of 5%. The amplitude of the furrows is 0.01 m (0.02 m peak-to-peak)

and their wavelength is 0.1 m. We choose a friction coefficient k = 0.04 m 1/3 s -1 . For the following computations, we use a time step ∆t = 0.001 s (this time step is imposed by the resolution of (2.1), as we need a small space step to get a good representation of the furrows).

Let us mention that all the numerical results are obtained using a C++ software for the resolution of the Shallow Water system, and a new library for the new friction coefficient.

Reference solutions and description of the test cases

This paragraph is devoted to the computation of the solutions of the Shallow Water model (2.1) where the geometry of the furrows is known explicitly. These solutions will be considered here as reference solutions. According to the parameters given above, the topography is modeled by the equation:

Z(x, y) = -0.05 y + 0.01 cos(20π y).

(4.10)

The space steps (with respect to x and y) are equal to 0.01 m, which means that each furrow is described by 200 cells. We assume that the domain is initially empty, that is u(0, x) = 0 and h(0, x) = 0. Let us denote by h, u, and q these reference solutions, at the small scale.

Rainfall test case

In this case, we impose rainfall on the whole domain with a constant permanent rain intensity R = 8 × 10 -4 m s -1 . The rain discharge is then Q R = 3.2 × 10 -3 m 2 s -1 . The final time is T = 22.5 s. Note that, since we are interested in the effects of the furrows, we focus on the transitional stage of the flow. Therefore the final time T is chosen such that the outflow discharge is approximately equal to the half of the rain discharge. We assume here that the upstream boundary is a solid wall. We show in Figure 5 the side-view of the water height at the final time.

Inflow test case

We also consider a permanent inflow from upstream. We prescribe

Q I = 3.132×10 -2 m 2 s -1
as discharge on the inflow boundary. The final time is T = 27.75 s. As for the rainfall test case, the final time was chosen such that the outflow discharge at T is approximately equal to the half of the inflow discharge. We show in Figure 6 the side-view of the water height at the final time. The space step with respect to y is set equal to the wavelength of the furrows, that is 0.1 m. The initial conditions remain unchanged, i.e. u(0, x) = 0 and h(0, x) = 0. We denote by capital letters (H, U and Q) the solutions of (2.4), at the large scale.

To improve the comparison in the rainfall test, we also compute the solution (h,u,q)

of the two-layer like system (2.5). The discretization parameters and the topography we use are the same as for system (2.4).

In order to compare the three models (2.1), (2.4) and (2.5), we consider the water height (h, H and h respectively) and the discharge (q, Q and q respectively) at the outflow. For this purpose, we first introduce h n i the average of the reference water height h contained in the furrow i at the time t n . We also consider H n i the water height in the furrow i at time t n computed with the model (2.4) for given K 0 and C, and H 0 n i in the case K 0 = 0 (that is the Shallow Water system on the coarser grid with plane topography and without the new friction term). Next we denote by e H the relative water height error defined by

e H =        N n=1 i h n i -H n i 2 N n=1 i h n i -H 0 n i 2        1/2 , ( 4.12)
where t N = T is the final time of the simulations. This error represents the effect of the new friction term K(h) in system (2.4) on the water height.

Finally we also study the discharge at the outflow: if q n is the value of q (from the resolution of (2.1) with the explicit topography) at time t n , for y = L, Q n and Q n 0 the solutions of (2.4) for given K 0 and C, and K 0 = 0 respectively, then the discharge error e Q is defined by:

e Q = N n=1 |q n -Q n | N n=1 |q n -Q n 0 | . (4.13)
The value of e Q shows the influence of the term K(h) on the discharge at the outflow.

Rainfall test case

For the rainfall test case (see paragraph 4.1.1), we present, in Figure 7, the water height error e H between models (2.1) and (2.4) as a function of K 0 for C = 10. We started with a step of 0.01 until K 0 = 0.1, and then the step is taken equal to 0.1 as the function increases.

We remark that the minimum of the error e H ≃ 0.2518 is obtained for K 0 = 0.02. Then by optimizing this error with respect to the two parameters K 0 and C, we finally find the minimum e H ≃ 0.1417 corresponding to K 0 = 0.02 and C = 0.4. The corresponding discharge error is e Q ≃ 5.8 × 10 -2 . We notice that the new model (2.4) allows to diminish the L 2 error on the water height by a factor 7 with respect to the case K 0 = 0, showing that the furrow effects are well taken into account.

We now report the results obtained with the two-layer like model (2.5), for the rainfall test. For different Manning's coefficients k, Figure 8 shows the water height error e h (see equality (4.12) replacing the solution H of system (2.4) by the solution h of (2.5)). The objective of this test is to understand how the changes in Manning's coefficient combined with a delay at the beginning of the experiment could represent the effects of the furrows.

We note that the minimum is e h ≃ 0.0422 for k = 0.03. The corresponding discharge error is e Q ≃ 8.6 × 10 -3 . We notice that the two-layer model allows to diminish the L 2 error on the water height by a factor 23 with respect to the model (2.4) for the case K 0 = 0.

The first conclusion on this rainfall test is that the two proposed models (2.4) and (2.5) are good representations of the effects of the furrows, as they lower the error to the reference solution compared to a plane topography. More precisely, the two-layer like model seems better on this test, but, as we mentioned before, its extension to the general two dimensional case is not obvious. To close this test case, we directly compare the results of the three models. Figure 9 gives the ratio between the outflow discharge obtained by the models and the rain discharge we imposed, as a function of time. We plot four results in the transitory regime: the discharge computed with (2.4) with K 0 = 0.02 and C = 0.4, the discharge with K 0 = 0, the solution of (2.5) with k = 0.03 and the reference solution q (see paragraph 4.1).

Although the domain in the new model (2.4) is an inclined plane, we observe that the additional friction term is able to retain the water for a moment during the rainfall. It is a good approximation to simulate the existence of the furrows. discharge ratio:

Q model /Q rain Shallow Water model new model with K 0 = 0 new model with K 0 = 0.02, C = 0.4 two-layer model simulated time (s)
Figure 9: Ratio between the outflow discharge and the rain discharge for all the models, as a function of time

Inflow test case

Now we consider the inflow test case (see paragraph 4.1.2) and we conduct the same study as for the rainfall test, in order to validate our new model of friction coefficient. Figure 10 shows the water height error e H between models (2.1) and (2.4) as a function of K 0 for C = 10. The step on K 0 is set equal to 0.001 when the error is decreasing, and 0.01 for larger errors. We notice that the minimum of the error is e H ≃ 0.3211 for K 0 = 0.004. We also present in Figure 11 the ratio between the outflow discharge Q n , Q n 0 and the imposed inflow at the transitory regime, for the model (2.4) with K 0 = 0.004, C = 10 and K 0 = 0 respectively and for the reference solution q computed in Section 4.1.

Here again we observe that the new model simulates well the effect of the furrows. Remark 4.1. In this part, we did not mention the optimization of the value of C: as the optimized coefficient for K 0 is ten times smaller than for the rainfall test, the results do not visibly vary when we change C.

Performance of the new model with calibrated coefficients

In the previous section, the numerical experiments show that our new model (2.4) with the additional friction coefficient to represent the furrows gives good results on the water height and the outflow discharge, compared to the reference solution computed on the explicit topography. But this model needs a calibration of the two coefficients K 0 and C.

In this paragraph, we study the robustness of this calibration when we change numerical (the space step) or physical (the slope or the roughness of the domain) parameters. For a future extension to fully two-dimensional problems, the new model has to predict the flow directions on each grid cell over the area of agricultural fields. As a consequence, the dimension of a grid cell can vary from 1 to 100 m 2 . So, a grid cell typically contains several furrows. To test for this scale effect, we present in Table 1 the relative errors e H and e Q to the reference solution as functions of the space step. This space step varies from 0.1 m to 0.4 m which means that a cell comprises 1 to 4 periods of the furrows. We can notice that the variation of the error e H is small: when the space step is multiplied by 4, the relative error is multiplied by less than 1.5.

Space step [m]

e At the same time, we present in Table 1 the computational times to compute the solutions of the new model (depending on the space step). The values are normalized with respect to the computational time to compute the reference solutions. As we can observe the new model allows to reduce the computational time for at least 90% and can easily achieve a 97% decrease under the condition that we tolerate a small error on the solution.

In Table 2, we present the errors e H and e Q for different general slopes of the topography (with a realistic range for agricultural fields), for the calibrated values K 0 = 0.02 and C = 0.4. We note that the new model allows to reduce the L 2 error on the water height by at least a factor 4 compared to the reference solution, showing that the effects of the furrows are well taken into account.

We complete these results with Table 3, where, for each slope, we calibrate the coefficients K 0 and C to get the minimum of the relative error e H . Comparing the values of the coefficients and of the errors in Tables 2 and3 .

Note that e H s has the same definition as e H (see (4.12)) except that, instead of summing the values of the water heights for all the discrete times, we only consider the time corresponding to the steady state t s .

The values we obtain with the new model are e H s ≃ 0.1765. Thus we deduce that at the steady state, the new model (2.4) with the calibrated coefficients takes well into account the effects of the furrows and approximates well the Shallow Water model (2.1).

Conclusions

In this paper, we proposed a new Shallow Water model (2.4) in order to describe the effects of furrows during overland flow without representing them explicitly. The main idea is to include in the classical Shallow Water equations (2.1) the additional friction term (2.3) that takes into account the effects of these furrows. The new model is proposed under the assumptions that the flow is fluvial and the direction of the flow is perpendicular to the furrows, parallel to a side of the rectangular domain. We also assumed that there is no infiltration and no soil erosion.

We presented numerical results to show the efficiency and the performance of the new model. We compared it with a two-layer like model, where there is a delay due to the filling of the furrows, and where we are able to assess the behavior of a model with a variable Manning's coefficient. Both models give good results for the rainfall test. We also showed that the calibration of the coefficient of the new model does almost not depend on the slope of the domain (see Table 2 and Table 3) and nor on the soil friction coefficient (see Table 4 and Table 5).

The numerical results presented in this paper are encouraging and indicate that the idea could be extended to more complex two-dimensional flows. This is now the main goal of the forthcoming works. Note that, unlike the two-layer model (2.5), the new model (2.4) can easily be generalized to two-dimensional problems. These extensions include the random variations of the height of the furrows. As we already mentioned in the description of (2.3), these variations can be taken into account in the constant C of the additional friction term. The extensions also include the case where the direction of the furrows is not perpendicular to the slope of the domain, in order to study the effects of this direction with respect to the slope.

Figure 1 :Figure 2

 12 Figure 1: Notations for a 1D Shallow Water flow

Figure 2 :

 2 Figure 2: An example of a topography with furrows that we consider

Figure 3 :

 3 Figure 3: Water trapped in a furrow

Figure 4 :

 4 Figure 4: Shape of the friction term K(h) for h F = 0.01 m

Figure 5 :

 5 Figure 5: Side-view of the water height at the final time for the rainfall test case

Figure 6 :

 6 Figure 6: Side-view of the water height at the final time for the inflow test case

Figure 7 :

 7 Figure 7: Water height error e H for the model (2.4) for the rainfall test case

Figure 8 :

 8 Figure 8: Water height error e h for the two-layer like model (2.5)

  new model with K 0 = 0 new model with K 0 = 0.004, C = 10 simulated time (s)

Figure 11 :

 11 Figure 11: Ratio between the outflow discharge and the imposed inflow discharge for all the models as a function of time

e

  Q for different Manning's coefficients T = t s = 50 s. The two models we consider are, in the one hand, the new model (2.4) with the calibrated coefficients K 0 = 0.02 and C = 0.4, and, on the other hand, the usual Shallow Water model (2.1) with an explicit description of the topography and of the furrows.We now denote e H s the water height error at steady state defined by

  The corresponding discharge error is e Q ≃ 4.4668 × 10 -2 . Note that our new model (2.4) allows to lower the L 2 error on the water height by a factor 1.62 with respect to the case K 0 .

	1.6
	1.4
	1.2
	1
	e H
	0.8
	0.6
	0.4
	0.2
	0	0.01	0.02	K 0	0.03	0.04	0.05
	Figure 10: Water height error e H for the model (2.4) for the inflow test case

Table 1 :

 1 Relative errors on the water height and the flowrate, and CPU times (normalized with respect to the time to compute the reference solution).

		H	e Q	CPU time
	0.1 0.2 0.4	0.1417 0.1675 7.2865 × 10 -2 5.8 × 10 -2 0.2024 8.963 × 10 -2	0.11 0.0574 0.0255

Table 2 :

 2 , we see that the new model with the calibrated values K 0 = 0.02, C = 0.4 still yields accurate results. Errors e H and e Q for different general slopes of the topography with K 0 = 0.02

	slope K 0	C	e H	e Q
	2%	0.02 0.4 0.2434	0.1855
	5% 8%	0.02 0.4 0.1417 5.8 × 10 -2 0.02 0.4 0.2194 0.187
	11% 0.02 0.4 0.2035	0.171
	and C = 0.4
	slope K 0	C	e H	e Q
	2%	0.02 0.3 0.2167	0.1134
	5% 8% 11% 0.04 0.4 0.1089 6.8554 × 10 -2 0.02 0.4 0.1417 5.8 × 10 -2 0.04 0.4 0.1205 8.1222 × 10 -2

Table 3 :

 3 Minimum errors e H with the corresponding coefficients K 0 , C and the values of e Q for different general slopes of the domain Table4shows the errors e H and e Q for different Manning's coefficients (the slope is kept equal to 5%), for the calibrated values K 0 = 0.02 and C = 0.4. As for the variation of the slope, in Table5, we give the coefficients K 0 and C that minimize the relative error e H for each Manning's coefficient, and we complete the table with the corresponding error e Q . Here again, we see that the results given by our model (2.4) are good compared to the reference solution.

	Manning's coefficient [m 1/3 s -1 ] K 0	C	e H	e Q
	0.001 0.04 0.1	0.02 0.4 0.0744 0.02 0.4 0.1417 0.02 0.4 0.2746 9.9691 × 10 -2 4.67 × 10 -2 5.8 × 10 -2

The last results we present are dedicated to water heights and discharges errors at steady state for the Shallow Water model (2.1). The final time of the simulations is then

Table 4 :

 4 Errors e H and e Q for K 0 = 0.02 and C = 0.4 for different Manning's coefficients

	Manning's coefficient [m 1/3 s -1 ] K 0	C	e H	e Q
	0.001 0.04 0.1	0.02 0.4 0.0744 4.67 × 10 -2 0.02 0.4 0.1417 5.8 × 10 -2 0.02 0.5 0.265 0.1128

Table 5 :

 5 Minimum errors e H with the corresponding coefficients K 0 , C and the values of

Acknowledgements

This work was supported by the ANR project METHODE #ANR-07-BLAN-0232 (http://www.univ-orleans.fr/mapmo/methode/) and the Région Centre, France. The authors wish to thank the members of the project for fruitful discussions and for many helpful comments.

• Compute U n+2 i using the Heun method defined by (3.8).

We now describe the treatment of the boundary conditions at the inflow and the outflow.

We denote by a b subscript the values on the (fictive) boundary cell, and by the index "in" the values in the first cell, inside the domain. The normal n is equal to -1 on the left boundary (x = 0) and 1 on the right boundary (x = ℓ).

• A solid wall is modeled imposing u b = -u in , h b = h in , on the condition that the topography be extended horizontally on the fictive cells.

• At the inflow boundary, we impose the discharge q b satisfying nq b < 0. Since we only consider fluvial flows, the water height h b is computed using Riemann invariants (see for instance [START_REF] Bristeau | Boundary conditions for the shallow water equations solved by kinetics schemes[END_REF][START_REF] Lhomme | One-dimensional, two-dimensional and macroscopic approaches to urban flood modelling[END_REF]). More precisely, assume that c = √ gh. It is well known that for the Shallow Water system (3.6), the quantity u ∓ 2c is constant along the characteristic dx dt = u ∓ c. Thus we have

Multiplying (3.9) by h b , we obtain

Newton method is used to solve this last equation and to get h b .

• At the outflow boundary, we always impose the water height h b and, here again, we use Riemann invariants to compute the discharge. We get

and we can easily deduce q b = h b u b .

All these previous steps are quite usual for the resolution of the Shallow Water system (2.1). Note that, for our new 2D model (2.4), the additional friction term is treated in an explicit way.

Numerical results

In section 2, we introduced several models: the usual Shallow Water system (2.1), the Shallow Water system with an additional friction coefficient that represents the furrows when we consider a plane topography (2.4) and a two-layer like model (2.5).

In this section, we present several results obtained with these three models in order to show the capacity of our new model (2.4) to approximate the exact solution. Namely, we


C Bohning

Fetecau

J Zierep

R Bohning

C Fetecau

Numerical results for the conservation of energy of maxwell media for the rayleigh-stokes problem

Keywords: Maxwell fluid, Dissipation, Power of the shear stresses at the wall, boundary layer thickness, Rayleigh-Stokes problem

. We start from the Fourier sine transform. The numerical result is given and discussed for the velocity u, the power of the wall shear stresses L, the dissipation Φ and the boundary layer thickness δ. These new results are important for nature and technology.

Introduction

In recent years, advancements in technological applications have brought a wide range of rheological complex fluids that are characterized by diverse significant deviations from the simple Newtonian behaviour. The applications of non-Newtonian fluids include extrusion of polymer fluids, exotic lubricants, colloidal and suspension solutions, food stuffs, synthetic propellants, molten plastics and many others. The departure from Newtonian behaviour manifests itself in a variety of ways [START_REF] Rajagopal | Mechanics of non-Newtonian fluids[END_REF]: non-Newtonian viscosity (shear thinning or shear thickening), stress relaxation, normal stress differences etc. Due to the diversity of fluids in nature many models have been proposed to describe their behaviour. Among them, the models of rate type as well as those of differential type have received much attention.

The first viscoelastic rate type model, which is still used widely, is due to Maxwell [START_REF] Maxwell | On the dynamical theory of gases[END_REF]. While Maxwell did not develop his model for polymeric liquids, but instead for air, the methodology that he used has been generalized by Rajagopal and Srinivasa [START_REF] Rajagopal | A thermodynamical frame-work for rate type models[END_REF] to produce a plethora of rate type models. However, Maxwell recognized that the body has a means for storing energy and a means for dissipating energy, the storing of energy characterizing the fluid elastic response and the dissipation of energy characterizing its viscoelastic nature.

The motion of a viscous fluid caused by the impulsive motion of a plane wall is termed as Stokes or Rayleigh´s problem [START_REF] Schlichting | Boundary Layer Theory[END_REF]. Stokes [5] solved the problem for a viscous fluid past an infinite plate. Soundalgekar [START_REF] Soundalgekar | Stokes problem for elastico-viscous fluid[END_REF] extended it to a fluid of second grade using a perturbation technique. Teipel [START_REF] Teipel | The impulsive motion of a flat plate in a visco-elastic fluid[END_REF] showed that for such a fluid a strict similarity solution does not exist and provided a series solution. Puri [START_REF] Puri | Impulsive motion of a flat plate in a Rivlin-Ericksen fluid[END_REF] and Bandelli et al. [START_REF] Bandelli | On some unsteady motions of fluids of second grade[END_REF] studied the problem and found solutions which do not satisfy the initial condition. This is not surprising due to the incompatibility of the prescribed data. However, their solutions are physically interesting and correct from mathematical point of view. A new exact solution for the same problem has been later established in [START_REF] Fetecau | On a class of exact solutions of the equations of motion of a second grade fluid[END_REF]. This solution, obtained by means of the Fourier sine transform, satisfy both the governing equation and all imposed initial and boundary conditions. Its extension to rate type fluids has been realized in [START_REF] Fetecau | A new exact solution for the flow of a Maxwell fluid past an infinite plate[END_REF] and [START_REF] Fetecau | The first problem of Stokes for an Oldroyd-B fluid[END_REF]. The velocity field and the adequate shear stress obtained in [START_REF] Fetecau | A new exact solution for the flow of a Maxwell fluid past an infinite plate[END_REF] have been recently used by Zierep and Fetecau [13] to develop a complete energetic study for the Rayleigh-Stokes problem of a Maxwell fluid. ____ Corresponding authors, E-mail address: * rainerbohning@gmx.de, ** cfetecau@yahoo.de

A c c e p t e d m a n u s c r i p t 2

Our purpose in this paper is to print out some new and interesting aspects regarding the behaviour of Maxwell fluids. More exactly, it comes as a complement to our previous work [START_REF] Zierep | Energetic balance for the Rayleigh-Stokes problem of a Maxwell fluid[END_REF] in which the approximate expressions as well as the diagrams of the three important physical entities L, Φ and δ ( the power of the shear stress at the wall, dissipation and the boundary layer thickness) have been evaluated for small values of the relaxation time, namely for 1 t / λ << . Here, these entities and their graphical illustrations have been evaluated for arbitrary values of t / λ . The differences are characteristic and very important for applications.

Statement of the problem and solution for the velocity u

We continue our publication [START_REF] Zierep | Energetic balance for the Rayleigh-Stokes problem of a Maxwell fluid[END_REF] in respect to the numerical results of the analytical solutions. We start with the Differential equation [13, (6)

] 2 2 2 2 y u ν t u λ t u ∂ ∂ = ∂ ∂ + ∂ ∂ . ( 1
) ) t , y ( u u =
-velocity, y, t -variables of length and time, ν -kinematic viscosity, λrelaxation time for elasticity.

The initial -boundary conditions for the Rayleigh -Stokes problem are [13, (10

)] 0 ) 0 , ( = y u , 0 > y , U t u = ) , 0 ( , 0 > t (2) 0 ) , ( = ∞ t u , 0 > t .
The Fourier sine transformation leads to ( )

( ) , sin 2 1 sin 1 1 2 1 cos sin 2 1 sinh 1 1 2 1 cosh 2 1 ) , ( 1 2 2 2 1 0 2 2 2 2 ⎟ ⎟ ⎟ ⎠ ⎞ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅ - - + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅ - + ⎜ ⎜ ⎜ ⎝ ⎛ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅ - - + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅ - ⋅ - = ∫ ∫ ∞ - ξ d ξ ξ θ η θ ξ ξ θ ξ ξ d ξ ξ θ η θ ξ ξ θ ξ e π U θ η u θ (3) t ν y η 2 = ; λ t θ = .
We get the numerical results for this Maxwell solution by suitable forms of integration.

Solution of the power of the wall shear stress L(t,λ)

Of great practical interest is the wall shear stress ) , 0 ( θ τ w and the power of the wall shear

stress U θ τ L w l ) , 0 ( = [ 13
, ( 14), ( 15)]

⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - = ∫ ∫ ∞ νλ νλ w ξ d γθ γ ξ d βθ β θ π μ θ τ U 2 1 2 1 0 2 sin 1 2 sinh 1 2 exp 4 ) ( 1 (4) 2 4 1 νλξ β - = , 1 4 2 - = νλξ γ .
Important are the two limits:

∞ → = θ λ t (Newtonian case) t πν U μ τ w - = , 0 → = θ λ t (solid body) νλ U μ τ w - = .
Fig. 2 shows the behaviour of ) (θ L . Of interest are again the limits

∞ → = θ λ t : t πν U μ U τ L w 2 l l - = = , 0 → = θ λ t : νλ U μ L 2 ) 0 ( l - = .
These limits agree completely with the numerical dates (Fig. 2).

Solution of the dissipation Φ(t,λ) and the boundary layer thickness δ(t,λ)

The calculation of Φ leads to a complicate way. First we need ) , ( t y τ and y

t y u ∂ ∂ ) , ( [13, (8)].
Then we have to split the integrals in a suitable way to get the numerical results. Fig. 3 gives the dissipation as function of t with the parameter λ. The behaviour is similar to the decay of a potential vortex under the influence of viscosity.

The boundary layer thickness [13, (9)] can be calculated from

∫ ∞ = 0 ) , ( ) , 0 ( 1 dy t y u t u δ
. We use here the velocity u(y, t) from chapter 2. Fig. 4 shows a typical example. Generally we make the following remark. The last three figures 2, 3 and 4 differ from the corresponding figures 1, 2 and 3 of the preceding publication [START_REF] Zierep | Energetic balance for the Rayleigh-Stokes problem of a Maxwell fluid[END_REF]. The reason for that is the

General discussion of the results

Today the energy balance for the movement of non-Newtonian fluids plays an important role.

We solved the Rayleigh-Stokes problem for Maxwell media. Starting from the analytical Fourier sine solution we have for the velocity u(y, t), the power of the wall shear stress L(t, λ) [13, ( 14), (15)], the dissipation ) λ , t ( Φ and the boundary layer thickness δ(t, λ) [13, (9)] a more or less difficult numerical problem to solve our Fourier sine analytical solution. This expensive method leads to the results given in Figs. 1234

Fig. 1 A c c e p t e d m a n u s c r i p t 3 some

 13 Fig. 1 shows a typical result. U u is given as function of t ν y η 2 = with λ t θ = as parameter.

A c c e p t e d m a n u s c r i p t 4 following

 4 . In[START_REF] Zierep | Energetic balance for the Rayleigh-Stokes problem of a Maxwell fluid[END_REF] the figures 1, 2 and 3 for L, Φ and δ have been calculated under the assumption 1 << t λ while in the above publication we calculate L, Φ and δ for arbitrary values of t λ . The differences are characteristic and important for applications.

  discussed above. The behaviour is interesting and leads immediately to the energy balance.

  Fig.1. U u as function of

Fig. 2 .

 2 Fig. 2. Power of wall shear stress L(t,λ) Glycerin: ν = 0.0011746 m 2 /s, μ = 1.48 s Pa , U = 0.5 m/s, l = 0.5 m.

A c c e p t e d m a n u s c r i p t 6 Fig. 3 .

 63 Fig. 3. Dissipation Φ(t) Glycerin: ν = 0.0011746 m 2 /s, μ = 1.48 s Pa , U = 0.5 m/s, l = 0.5 m.

Fig. 4 .

 4 Fig.4. Boundary layer thickness δ(t,λ) [m],Glycerin: ν = 0.0011746 m 2 /s.


T Károlyi

Bódai

A Bibó

Gy Károlyi

T Bódai

Keywords: Bouncing ball, Elastic impact, Oblique impact

The to and fro motion of a bouncing ball on a flat surface is represented by a low-dimensional model. To describe the repeated reversals of the horizontal velocity of the ball, the elasticity of the ball has to be taken into account. We show that a simple fly-wheel model exhibits the observed hither and thither motion of elastic balls. The suggested model is capable of describing oblique impacts of spherical bodies, which can be important in many applications, including dynamical simulation of granular materials. We find that the behaviour of the bouncing fly-wheel is sensitive to the initial conditions, and the escape time plots are used to illustrate this observation.

A c c e p t e d m a n u s c r i p t

Fly-wheel model exhibits the hither and thither motion of a bouncing ball

Introduction

Collisions are very important in many fields of Nature and Science and have wide-spread applications in industry. In Nature, the impact of ice particles are thought to be related to the formation of Saturn's rings [START_REF] Bridges | Structure, stability and evolution of Saturn's rings[END_REF], and the singleparticle impacts can play a crucial role in understanding the behaviour of snow [START_REF] Schaer | Particle densities, velocities and size distributions in large avalanches from impact-sensor measurements[END_REF][START_REF] Bartelt | Avalanche physics ploughs ahead[END_REF] or sand [START_REF] Andreotti | The song of dunes as a wave-particle mode locking[END_REF] avalanches. In Science, the ergodic behaviour of impacting billiard systems is supposed to give ideas to describe the stochasticity of gas and to the Boltzmann hypothesis [START_REF] Dorfman | An introduction to chaos in nonequilibrium statistical mechanics[END_REF]. In many sports like tennis, golf or baseball it takes years of practice to achieve the desired bounce of the ball with the right speed, angle and spin [START_REF] Cross | The coefficient of restitution for collisions of happy balls, unhappy balls, and tennis balls[END_REF][START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF][START_REF] Cross | Experimental study of the gear effect in ball collisions[END_REF]. Applications of impacting bodies include the use of granular materials in chemical, pharmaceutical, agricultural, mining

A c c e p t e d m a n u s c r i p t

or mineral processing industry [START_REF] Dong | Measurement of impact behaviour between balls and walls in grinding mills[END_REF][START_REF] Dong | Experimental study of oblique impacts with initial spin[END_REF]. In computer simulations of granular materials [START_REF] Sondergaard | Measurement of solid speheres bouncing off flat plates[END_REF][START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Brilliantov | Model for collisions in granular gases[END_REF][START_REF] Labous | Measurements of collisional properties of spheres using high-speed video analysis[END_REF][START_REF] Louge | Anomalous behavior of normal kinematic restitution in the oblique impacts of a hard sphere on an elastic plate[END_REF][START_REF] Vu-Quoc | An accurate tangential forcedisplacement model for granular-flow simulations: Contacting spheres with plastic deformations, force-driven formulation[END_REF] it is desirable to have simple, but accurate models of the interaction forces and impact properties.

In order to construct such models many experimental investigations of impacts have been carried out. Probably the simplest of these experiments concern the oblique collision of balls or disks with a planar surface, with or without initial spin [START_REF] Cross | The coefficient of restitution for collisions of happy balls, unhappy balls, and tennis balls[END_REF][START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF][START_REF] Cross | Experimental study of the gear effect in ball collisions[END_REF][START_REF] Dong | Measurement of impact behaviour between balls and walls in grinding mills[END_REF][START_REF] Dong | Experimental study of oblique impacts with initial spin[END_REF][START_REF] Sondergaard | Measurement of solid speheres bouncing off flat plates[END_REF][START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Labous | Measurements of collisional properties of spheres using high-speed video analysis[END_REF][START_REF] Louge | Anomalous behavior of normal kinematic restitution in the oblique impacts of a hard sphere on an elastic plate[END_REF][START_REF] Maw | The rebound of elastic bodies in oblique impact[END_REF][START_REF] Maw | The role of elastic tangential compliance in oblique impact[END_REF]. Parallel with the experimental studies, many models have been formulated to describe oblique impacts. The early attempts considered the impacting ball to be rigid. In Brody's model [START_REF] Brody | That's how the ball bounces[END_REF] the ball either slides along the surface throughout the duration of the impact, or it stops sliding due to the friction force and rolls without slip thereafter. However, such rigid body models are not capable of describing the observed changes in the direction of the tangential force acting on the contact surface [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF][START_REF] Cross | Experimental study of the gear effect in ball collisions[END_REF] or reversals of the contact point velocity [START_REF] Dong | Measurement of impact behaviour between balls and walls in grinding mills[END_REF] during a single impact of a ball, observed most notably in case of solid rubber balls referred to as "superballs" [START_REF] Johnson | The bounce of 'superball[END_REF][START_REF] Stronge | Impact Mechanics[END_REF][START_REF] Garwin | Kinematics of an ultraelastic rough ball[END_REF]. In Brody's model once rolling starts, no mechanism can re-initiate sliding during the single collision. To avoid this problem, several models have been developed, based either on the assumption of a no-slip contact point and an elastic tangential force [START_REF] Garwin | Kinematics of an ultraelastic rough ball[END_REF], or on a prescribed number of reversals in the slip direction of the rigid ball [START_REF] Keller | Impact with friction[END_REF][START_REF] Stronge | Rigid body collisions with friction[END_REF].

It became clear that the elastic deformations must be taken into account for a proper modelling of the tangential forces. In principle, the deformable surface of the impacting ball acts like a "tangential spring" during the collision with the surface [START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Maw | The rebound of elastic bodies in oblique impact[END_REF][START_REF] Maw | The role of elastic tangential compliance in oblique impact[END_REF][START_REF] Maw | The oblique impact of elastic spheres[END_REF]. The elastic deformations of balls compressed to each other in the presence of a tangential force and a torsional couple were first described by Mindlin [START_REF] Mindlin | Compliance of elastic bodies in contact[END_REF], Mindlin and Deresiewicz [START_REF] Mindlin | Elastic spheres in contact under varying oblique forces[END_REF], and Walton [START_REF] Walton | The oblique compression of two elastic spheres[END_REF]. They have found that the deformations and the force distribution on the contact surface depend on the complete quasi-static loading history. They have revealed that the contact surface is divided into annuli that are either in the state of slip or no-slip, the slip regions being on the outer annulus of the contact surface. Based on these essentially static considerations, Maw et al. have constructed a model [START_REF] Maw | The oblique impact of elastic spheres[END_REF] of collision where the stick and slip regions coexist at each time instant during the impact. In this model, the contact surface consists of the same slip and no-slip annuli as in Mindlin's static model, but they change with time. This makes it necessary to find by an iteration in each time-step which of the annuli are in slip or in no-slip, which makes the solution a cumbersome process. However, the predictions of the model coincide remarkably well with experimental observations [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF][START_REF] Dong | Experimental study of oblique impacts with initial spin[END_REF][START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Louge | Anomalous behavior of normal kinematic restitution in the oblique impacts of a hard sphere on an elastic plate[END_REF][START_REF] Maw | The rebound of elastic bodies in oblique impact[END_REF][START_REF] Maw | The role of elastic tangential compliance in oblique impact[END_REF].

It is the interplay of the repeated stick and slip [START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Maw | The oblique impact of elastic spheres[END_REF] and the effects of the elastic deformations [START_REF] Johnson | The bounce of 'superball[END_REF][START_REF] Maw | The oblique impact of elastic spheres[END_REF] that causes the complexity of the oblique impact of a ball on a flat surface. As a consequence of this, with some practising, one can throw a ball in such a way, that it "comes back" after the collision with the surface. The ball has to spin "backwards", as shown in Fig. 1, to achieve this reversal. In fact, it is possible to throw a solid rubber ball (e.g. superball) so that in the first few bounces on a flat horizontal surface it always rebounces, that is, the velocity of its centre of mass and its angular velocity changes direction in

A c c e p t e d m a n u s c r i p t

each collisions [START_REF] Dong | Experimental study of oblique impacts with initial spin[END_REF][START_REF] Maw | The role of elastic tangential compliance in oblique impact[END_REF][START_REF] Johnson | The bounce of 'superball[END_REF][START_REF] Stronge | Rigid body collisions with friction[END_REF]. In this paper we address this hither and thither motion of an elastic ball, and construct a very simple, low-dimensional model that can exhibit a very similar behaviour.

v ω Figure 1: Collision of the ball with the ground with back-spin. The horizontal velocity of the lowest point of the ball is larger than v, that of the centre of mass, due to the direction of the angular velocity ω. The friction slows the lowest point until it stops, then rolling would set in for a rigid ball.

The only force that can explain why the ball changes its horizontal velocity at the collision is the friction force. To change the momentum of the ball, the friction force must act for a finite, nonzero period of time. This collision time will be denoted by τ . During this time the friction force F (t) exerts a time-dependent torque on the ball, and hence changes its angular momentum. Imagine that the ball collides with the ground while spinning "backwards" as shown in Fig. 1, so that its lowest point moves faster than its centre of mass. To achieve the to and fro motion, the ball has to change the direction of both its velocity and angular velocity during the collision. However, in order to move and rotate in the opposite direction, the contact surface (or contact point, in case of a rigid ball) has to stop for an instant, which means that the sliding motion during the collision ceases. For a rigid ball, in the absence of other effects, the friction force between the ball and the surface disappears, and from then on the ball would never slide again, its rotation and translation would synchronise. The rigid ball cannot change the directions of both its velocity and angular velocity. Hence in order to explain the back and forth motion of the ball, it cannot be considered rigid, deformations has to be incorporated in the model. We construct a low dimensional model that can give account of the role of deformations, and can exhibit the repeated rebounces.

In Section 2 we introduce the fly-wheel model that is capable of modelling deformations by adding a new degree of freedom to the rigid ball, characterising the internal torsional deformations. Then in Section 3 we analyse the forces acting during the impact, following the theory of Hertz [START_REF] Hertz | Über die Berührung fester elastische Körper (On the contact of elastic solids)[END_REF], and compare the behaviour of the fly-wheel model with previous experimental findings. In Section 4 the equations of motion are cast into a dimensionless form, then Section 5

A c c e p t e d m a n u s c r i p t

contains numerical simulations of the motion of the fly-wheel model exhibiting the hither and thither motion. Then in Section 6 we summarise our results.

The fly-wheel model of the bouncing ball

The fly-wheel model incorporates the global effects of the internal elastic torsion of the ball in the simplest possible way. To account for the difference in the deformation of the outer part slowed down by friction and that of the inner part, we construct a simple fly-wheel model (Fig. 2). It consists of two parts, for simplicity shown as disk-like objects, that are allowed to rotate around a common axis, while they are connected by elastic torsional springs. The outer cylinder with radius R o mimics the surface layer of the ball, while the inner one with radius R i models the inner parts of the ball. During a collision, similar to the one depicted in Fig. 1, the outer part might stop due to friction, but the inner part behaves as a fly-wheel storing energy. Through the torsional springs it can exert an internal torque M on the outer cylinder, which can start moving the outer part again. We show that the internal elastic rotations can repeatedly reverse the horizontal velocity of the fly-wheel with the help of the friction force. The position of the centre of mass of the fly-wheel, that is, the common axis of the two disks is described by coordinates x(t) and y(t). We assume that the

A c c e p t e d m a n u s c r i p t

movement takes place in the x, y plane, and the axis is perpendicular to this plane. The vertical coordinate y is positive upwards, and x is the horizontal axis. The origin is chosen such that y = 0 when the outer disk just reaches the surface at the beginning and at the end of the collisions. The rotation of the outer disk is ϕ o (t), that of the inner is ϕ i (t). The internal torque in the torsional spring is taken to be

M (t) = k(ϕ o (t) -ϕ i (t)),
where k is the spring constant. Now we describe the motion of the fly-wheel model from the end of one impact when free flight starts, to the next such instant. We introduce the notion () f for quantities just after a collision when free flight starts, and () c for quantities just before the next collision occurs.

We assume that the internal relative rotations are damped between two collisions. Before the next collision the angles and the angular velocities of the two disks become the same: ϕ oc = ϕ ic and φoc = φic , where the dot indicates derivation with respect to time. During the flight, however, the total angular momentum is preserved, which means that the common angular velocity of the disks must be φoc

= φic = Θ o φof + Θ i φif Θ o + Θ i (1)
before the following collision. Here Θ o and Θ i are the moments of inertia of the outer and inner disks, respectively. Assuming strong dissipation of the internal rotations during the free flight, we can take

ϕ oc = ϕ ic = Θ o ϕ of + Θ i ϕ if Θ o + Θ i + φoc t fl , (2)
where t fl = 2 ẏf /g is the duration of the free flight. Between two collisions the centre of mass of the fly-wheel moves under the action of the gravitational force.

The position and velocity components of the fly-wheel right after a flight phase can be computed from the same data given at the beginning of the flight:

x c = x f + 2 ẋf ẏf g , y c = y f = 0, ẋc = ẋf , (3) ẏc
= -ẏf , ϕ oc = ϕ ic = Θ o ϕ of + Θ i ϕ if Θ o + Θ i + 2 ẏf φoc g , φoc = φic = Θ o φof + Θ i φif Θ o + Θ i ,
where g = 9.81 m/s 2 is the gravitational acceleration. The collision takes place during a short time τ . During this time, the forces and torques acting on the inner and outer cylinder are shown in Fig. 3. Note that during the impacts the gravitational force is generally much smaller than the impact force, hence it is neglected [START_REF] Keller | Impact with friction[END_REF]. It is also assumed for simplicity

A c c e p t e d m a n u s c r i p t

that the contact surface is very small, point-like, and hence no couple occurs there. This assumption, strictly true only for low impact speed, simplifies the analysis, a treatment of the effects of a couple resulting from an extended contact surface can be found in Ref. [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF]. The horizontal acceleration of the fly-wheel is ẍ = F/m, where m is the total mass of the system, and F = F (t) is the frictional force. The force A acting between the inner and outer parts via their common axis is also shown, but, being an internal force, it does not appear in the equations of motion. The inner cylinder's rotation is governed by the moment M imposed by the torsional springs: Θ i φi = M . For the outer disk, we have

Θ o φo = -M + F R o .
The moment in the springs depends on the relative rotation of the disks:

M = k(ϕ o -ϕ i ),
where k is the spring constant.

Introducing φ = ϕ oϕ i as the angle difference between the outer and inner parts leads to the following system of equations:

ẍ = F (t)/m, φ ≡ φo -φi = - k Θ o + k Θ i φ + F (t)R o Θ o , φo = 1 Θ o [F (t)R o -kφ], φi = k Θ i φ. (4)
Note that the equations contain a non-smooth term, the F (t) frictional force, which is different in case of sliding (slipping) and rolling (sticking).

In order to compute the friction force F (t), we must know whether the system is sliding or rolling. In fact, during one single impact, there might be transitions between sliding and rolling [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF][START_REF] Dong | Experimental study of oblique impacts with initial spin[END_REF][START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Louge | Anomalous behavior of normal kinematic restitution in the oblique impacts of a hard sphere on an elastic plate[END_REF][START_REF] Maw | The rebound of elastic bodies in oblique impact[END_REF][START_REF] Maw | The role of elastic tangential compliance in oblique impact[END_REF]. When there is sliding, the friction force is F = -μN sgn( ẋ + R o φo ), where N is the vertical force between the ground and the outer disk, μ is the friction coefficient, and sgn(a) is the sign of its argument: it is 1 if a is positive, it is -1 for a < 0, and for a = 0 it is 0.

Rolling occurs when the contact point stops, ẋ + R o φo = 0. From the first derivative of this with respect to time, substituting ẍ, the friction force turns out to be F = -mR o φo = mR o kφ/(Θ o + mR 2 o ). In case of rolling the friction force is proportional to φ, and hence to the internal moment M . Contrary to rigid balls, because of moment M the friction force can be different from zero in case of rolling.

Summing up the equations valid during time τ of the bounce, we thus have two distinct cases. First, the equations for the case of sliding are:

ẍ = - 1 m μN sgn( ẋ + R o φo ), φ = -k 1 Θ o + 1 Θ i φ - R o Θ o μN sgn( ẋ + R o φo ), φo = - k Θ o φ - R o Θ o μN sgn( ẋ + R o φo ),
F . φi = k Θ i φ. (5)
The equations valid during rolling are:

ẍ = R o k Θ o + mR 2 o φ, φ = -k 1 Θ o + mR 2 o + 1 Θ i φ, φo = - k Θ o + mR 2 o φ, φi = k Θ i φ. (6)
Note that in both cases the dynamics of ϕ i , the inner rotation, just follows φ, while in case of rolling all variables follow the dynamics of φ. The transition from sliding to rolling occurs when the velocity of the contact point becomes zero, ẋ + R o φo = 0, [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF] the transition from rolling to sliding occurs when the force required to provide

A c c e p t e d m a n u s c r i p t

the horizontal acceleration reaches the maximum available friction force:

|F | = M mR o Θ o + mR 2 o ≤ μN. (8)
Here we assume that the coefficient of sliding and static frictions are the same.

For simplicity, the fly-wheel model consists of disks, but the model would essentially be unaltered for a ball as well, only the moments of inertia would need to be chosen accordingly. In order to solve Eqs. ( 5) and ( 6), we need to know how N , the vertical force between the ground and the outer cylinder depends on time. We also need to know the duration τ of the collision. To find these, we assume that the collision is elastic, and use the results of Hertz's theory [START_REF] Hertz | Über die Berührung fester elastische Körper (On the contact of elastic solids)[END_REF][START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF][START_REF] Landau | Theory of elasticity[END_REF]. This is detailed in the next section.

Collision

Probably the simplest theory to describe the normal component of the collision of two elastic bodies is that of Hertz [START_REF] Hertz | Über die Berührung fester elastische Körper (On the contact of elastic solids)[END_REF]. His theory is limited to cases when the elastic oscillations are negligible, and the impact velocity is smaller than the velocity of sound (low frequency vibrations are absent) [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF][START_REF] Landau | Theory of elasticity[END_REF][START_REF] Hunter | Energy absorbed by elastic waves during impact[END_REF]. Also, we assume that the deformations during the collision remain elastic [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF][START_REF] Landau | Theory of elasticity[END_REF], and the effect of external forces (like gravity) can be neglected compared to the impact force [START_REF] Keller | Impact with friction[END_REF]. The collision time proposed by this theory has been confirmed in experiments [START_REF] Bokor | The measurement of initial impact velocity and contact time[END_REF], and it seems to be in agreement with measurements well beyond the validity of the assumptions. Following Hertz's theory of collision, assume that the centres of two elastic balls of radius R and R are at a distance R + Rh(t) from each other during their collision. Here h(t) is the effect of the deformation, and is related to the force N (t) acting between the two balls as [START_REF] Hertz | Über die Berührung fester elastische Körper (On the contact of elastic solids)[END_REF][START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF][START_REF] Landau | Theory of elasticity[END_REF]

h(t) = N (t) 2/3 D 2/3 (1/R + 1/R ) 1/3 , (9)
where

D = 3 4 [(1 -σ 2 )/E + (1 -σ 2 )/E ],
E and E are Young's moduli, and σ and σ are Poisson's ratios of the materials of the two balls. Modelling the collision of the ball with a rigid plane as R → ∞, E → ∞, and identifying R with R = R o leads to

h(t) = N (t) 2/3 D 2/3 R -1/3 o with D = 3(1 -σ 2 ) 4E . ( 10
)
If the mass of the ball is m, and we introduce

κ = 4 √ R o 5D , (11)
the conservation of energy for the elastic collision can be written as [START_REF] Hertz | Über die Berührung fester elastische Körper (On the contact of elastic solids)[END_REF]

m ḣ2 + κh 5/2 = m ẏ2 c , (12)

A c c e p t e d m a n u s c r i p t

where ẏc is the velocity of the ball before the collision. With ḣ = 0 we obtain the maximum deformation

h max = m ẏ2 c κ 2/5 , (13)
and integrating [START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF] between 0 and h max gives τ/2, the half of the collision time:

τ 2 = hmax 0 dh ẏ2 c -κ m h 5/2 . ( 14
)
This leads to [START_REF] Hertz | Über die Berührung fester elastische Körper (On the contact of elastic solids)[END_REF][START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF][START_REF] Landau | Theory of elasticity[END_REF]]

τ = 2 √ π Γ(7/5) Γ(9/10) m 2 -ẏc κ 2 1/5 (15)
as the duration of the impact, where Γ() is the Gamma function.

The deformation h(t) time t after the start of the collision can also be computed by integrating [START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF] between deformations 0 and h, resulting in

t = h 0 dh ẏ2 c -κ m h 5/2 = τ 2 √ π Γ 9 10 Γ 7 5 h(t) h max 2 F 1 2 5 , 1 2 , 7 5 , h(t)
h max 5/2 , (16)
where 2 F 1 () is the (Gaussian) hypergeometric function [START_REF] Luke | The special functions and their approximations[END_REF][START_REF] Bell | Special functions for scientists and engineers[END_REF]. The graph of this function is shown in Fig. 4 with solid line.

The dependence of the deformation on time, as given by ( 16), can be very closely modelled by a sine function. To make computations easier, instead of Eq. ( 16) we use the approximation

h(t) = h max sin πt τ . (17)
This approximation, used already in Ref. [START_REF] Hunter | Energy absorbed by elastic waves during impact[END_REF], is illustrated in Fig. 4, with a very good agreement. We can thus use this expression to compute the vertical force between the outer disk and the ground during the collision with the aid of [START_REF] Dong | Experimental study of oblique impacts with initial spin[END_REF] as

N (t) = R o h(t) 3 D . ( 18
)
As we assume completely elastic collisions, the vertical velocity after the impact is simply the opposite of that before the collision: ẏf =ẏc .

(

) 19
We also assume that the tangential force component does not affect the normal component, therefore we use the Hertz model for the normal component of oblique impacts. 16). The time is measured in units of the total time τ of the collision, the deformation is measured relative to the maximum deformation h max . The approximation according to Eq. ( 17) is also shown with dotted line.

The fly-wheel model was tested with parameters set to mimic the bouncing "superball" used in experiments by Cross [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF]. In that paper, among other types of balls, a superball of radius R o = 2.3 cm and mass m = 46.4 g has been used, and its bouncing properties have been measured. We estimate the friction coefficient to be μ = 0.8 (typical value for rubber on e.g. concrete [36]). The Poisson ratio for rubber is very close to σ = 0.5 [START_REF] Rinde | Poisson's ratio for rigid plastic foams[END_REF], and the elastic modulus is approximated with E = 2 MPa [START_REF] Hill | Large elastic compression of finite rectangular blocks of rubber[END_REF]. As an estimate, we take Θ o = (Θ o + Θ i )/5 = 1 5 • 2 5 mR 2 o , which gives Θ o = 19.64 g•cm 2 and Θ i = 78.55 g•cm 2 with the mass and radius given above. Fixing the internal disk, a force F acting on the ball's surface causes a rotation φ in the fly-wheel model, hence k = F R o /φ holds. Assuming that the rotation is the consequence of shear deflection in case of the superball, we find that φ = F/2GR 2 o , where F/2R 2 o is the approximated shear stress from the force F . This implies that the spring constant is k = 2R 3 o G, where the shear modulus is G = E/3 for rubber with σ = 0.5. This simple estimate gives k = 2 3 R 3 0 E = 16.2 Nm. With these parameter values, Fig. 5 shows the normal and tangential forces during an impact computed from our fly-wheel model. The initial conditions (impact velocity of v = 2.69 m/s at an angle 36 • with the surface, without spin) approximate those in the experiment by Cross [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF]. The obtained forces are very close to those measured in Ref. [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF], even with our crude estimation of the parameters. Similar curves have been obtained for other types of balls as well [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF]. As a further test of our model, we show in Fig. 6 how the tangential velocity of the contact point changes in the impact. In Fig. 6 the tangential velocity of the contact point before the impact is described by the dimensionless β c = -( ẋc +R o φoc )/μ ẏc , the tangential velocity after the impact by β f = -( ẋf + R o φof )/μ ẏc . The curve, obtained for impacts without initial spin ( φoc = 0) follows closely the experimental results of Ref. [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF], and is similar to those obtained for different types of balls or colliding disks [START_REF] Cross | Grip-slip behavior of a bouncing ball[END_REF][START_REF] Dong | Experimental study of oblique impacts with initial spin[END_REF][START_REF] Schäfer | Force schemes in simulations of granular materials[END_REF][START_REF] Louge | Anomalous behavior of normal kinematic restitution in the oblique impacts of a hard sphere on an elastic plate[END_REF][START_REF] Maw | The rebound of elastic bodies in oblique impact[END_REF][START_REF] Maw | The role of elastic tangential compliance in oblique impact[END_REF][START_REF] Johnson | The bounce of 'superball[END_REF].

We note that Stronge et al. have developed a similar low-dimensional model for oblique impacts [START_REF] Stronge | Oblique impact with friction and tangential compliance[END_REF]. In their model, a massless "particle", connected by horizontal and vertical springs to the rigid bouncing body, makes the connection with the surface. This particle with the springs models the compliance of the impacting region. Due to the zero mass of the contacting particle, in their model they need to specify whether there is initial sliding or rolling, as the massless particle can stop immediately. This is expected to be similar to the behaviour of our fly-wheel model for Θ o = 0, that is, when the inertia of the outer disk is zero. They have also found the reversal of the tangential force in their model. The dimensionless tangential velocity β f of the contact point after the impact as a function of that before the impact (β c ). The parameters used for the simulation are the same as those used for Fig. 5, but the velocity before the impact was always 4 m/s. For small incident angles (measured between the direction of the incoming velocity and the surface normal, β c 7), the ball starts to roll during the impact, for large angles (β c 7) it slides throughout the contact.

Dimensionless equations

The governing equations derived so far can be cast into a non-dimensional form. Let

= Θ o /Θ i , ξ = 2k/R o mg, ζ = Θ o /mR 2 o and α = kD/R 3 o .
Omitting the tildes from the new dimensionless variables for brevity, we can rewrite Eq. ( 3) describing the free flight as 2/5 α 2/5

x c = x f + ξ ẋf ẏf , y c = y f = 0, ẋc = ẋf , (20) ẏc
= -ẏf , ϕ oc = ϕ ic = νϕ of + ϕ if ν + 1 + ξ ẏf φoc , φoc = φic = ν φof + φif ν + 1 , φ c = ϕ oc -ϕ ic , φc = φoc -φic ,
ζ 3/5 (-ẏc ) -1/5 , (22)
while the maximum deformation (13) becomes

h max = 5 4 α 2/5 ζ -3/5 ẏ4/5 c . (23)
The dimensionless vertical force [START_REF] Maw | The role of elastic tangential compliance in oblique impact[END_REF] during the collision, using [START_REF] Maw | The rebound of elastic bodies in oblique impact[END_REF], can be computed as

N (t) = 1 α h max ζ sin πt τ 3/2 , (24)
where time t is measured from the start of the collision. This force is used to compute the course of the collision. When the friction force F exceeds μN , that is, when

|φ| > (1 + ζ)μN, (25)
However, when the velocity of the lowest point, in contact with the ground, becomes zero, that is,

ζ ẋ + φo = 0, (27)
rolling starts, and the dimensionless form of the equations of motion (6) become

ẍ = 1 1 + ζ φ, φ = - ζ 1 + ζ + ν φ, φo = - ζ 1 + ζ φ, φi = νφ. (28)
After the collision, the vertical velocity changes sign:

ẏf = -ẏc . (29)

A c c e p t e d m a n u s c r i p t

The dimensionless parameters α, ζ, ξ, ν and μ characterise the system. However, if we are to model a bouncing ball, not all of these parameters are independent. Because the total moment of inertia of a ball is Θ o + Θ i = 2 5 mR 2 0 , we find the relation

ζ = 2 5 • ν 1 + ν (30)
between the dimensionless parameters. (For a disk-like fly-wheel with Θ o +Θ i =

which implies ζ = 0.08, according to [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF]. These values correspond to the parameters used in the previous section to describe the bounce of the superball, except that we use k = 4.0 Nm for the spring constant.

Simulation results

Now we have all the equations to compute the movement of the bouncing flywheel. We use Eqs. [START_REF] Johnson | The bounce of 'superball[END_REF] to to find the state of the system before the collision. Then Eq. ( 24), together with [START_REF] Keller | Impact with friction[END_REF], gives the vertical force N (t) as a function of time during the collision, and Eq. ( 22) gives the duration τ of the impact. Knowing these, we can solve numerically the dynamical equations ( 26) or ( 28), depending on whether the contact between the outer disk and the ground is sliding or rolling, respectively. The new vertical velocity after the impact comes from [START_REF] Hertz | Über die Berührung fester elastische Körper (On the contact of elastic solids)[END_REF]. To follow the motion numerically, the 4th order Runge-Kutta method with a fixed time-step of τ/10 5 is applied during the impacts of length τ . This resolution proves to be sufficient, further refinement does not change the trajectories significantly. The initial values for the coordinates x 0 and y 0 of the ball, its initial velocity components v x0 , v y0 , and the initial spin ω 0 ≡ φo (0) = φi (0) of the fly-wheel need to be set (we assume that at the release of the ball the inner and outer parts spin together). We use an initial dimensionless release height of y 0 = 543.48 (corresponding to 1 m) and initial horizontal position x 0 = 0. We assume that at release the ball has only horizontal velocity, v y0 = 0; this is not an important condition, to achieve different values of the vertical velocity we could release the ball from a different height.

In Fig. 7 an example trajectory is shown, where the fly-wheel performs the hither and thither motion. The fly-wheel was released with a considerable backspin, and the ball was bouncing back and forth for a long time.

In order to get a global picture of the fly-wheel's behaviour, the dependence on the initial conditions is surveyed in Fig. 8. Here, a colour coding is used to indicate how many times the fly-wheel, released with the corresponding initial conditions, changes its horizontal velocity in the impacts before settling, that is, before bouncing away in the same direction. The lighter blue colours indicate many rebounces, the lighter red areas correspond to a small number of rebounces. We see that the initial conditions leading to a large number of rebounces accumulate along a straight line. As we approach this line, the number of rebounces increases rapidly. The equation of the line, for this parameter setting, is v x0 = 5ω 0 . In Fig. 8 the parity of the number of reversals changes when we cross the line of high number of rebounces. This means that the line separates regions where the direction of the final bouncing away is different, and at this line the fly-wheel is "hesitating" which way to go after the hither and thither motion. In this respect the line separates the initial conditions leading to qualitatively different outcomes, and the process is very similar to scattering processes [START_REF] Tél | Chaotic dynamics: An introduction based on classical mechanics[END_REF][START_REF] Bleher | Routes to chaotic scattering[END_REF][START_REF] Ott | Chaotic scattering: An introduction[END_REF][START_REF] Jung | Scattering one step from chaos[END_REF].

To test the effect of the ratio ν of the moment of inertia of the inner to that of the outer part, we show in Fig. 9 how the number of reversals depends on ν and the initial spin ω 0 . Note that, according to [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF], if the value of ν is varied, the value of ζ changes as well. The extremal value ν = 0 (and accordingly, ζ = 0) denotes the case when the outer part has no moment of inertia, and hence the fly-wheel consists only of the internal disk with a massless outer part (very strong fly-wheel effect). In this range we find high number of rebounces with a small initial spin. For large ν, that is, when the moment of inertia of the outer disk is large, the energy stored in the internal part can be insufficient to turn back the rotation of the outer disk for small initial spin, the necessary backspin to see the hither and thither motion increases. [START_REF] Landau | Theory of elasticity[END_REF], further initial conditions are x 0 = 0, ẏ(0) = 0, y 0 = 543.48, and 500 × 500 initial conditions have been followed for a).

Summary

"Angle of incidence equals angle of reflection" is a common concept for bouncing balls. However, throwing a real ball with backspin, it can reverse its horizontal velocity several times before it bounces off in one direction. Even if considered in two dimensions, besides translation the ball might possess a spinning motion. This implies that the point of contact may move at a speed different from that of the centre of gravity, and it leads to a friction force acting between the ball and the surface that can turn it back during the impact of finite (non-zero) duration. This effect alone is not enough, however, to explain the repeated reversals in the horizontal velocity, the elasticity of the ball has to be taken into account as well.

A simple model, the fly-wheel, takes into account both of these important effects. The stick and slip intervals during impacts provide the necessary friction force, while the internal rotation of the fly-wheel models the tangential elastic deformations throughout the collisions. We have shown that this simple model is capable of providing realistic impact and friction forces during collisions, and the relationship between incident and reflection velocities is close to that experimentally observed. This model, because of its simplicity, might prove to be useful for simulations of granular materials or other similar applications.

The fly-wheel model can exhibit the hither and thither motion of bouncing balls. This motion, illustrated in Fig. 7 is similar to that observed in wedge billiards [START_REF] Lehtihet | Numerical study of a billiard in a gravitational field[END_REF][START_REF] Tél | Chaotic dynamics: An introduction based on classical mechanics[END_REF]. In that case, an ideal, point-like ball is bouncing between two intersecting inclined planes, where the impacts are assumed to follow the simple reflection rule: angle of incidence equals angle of reflection. The wedge billiard exhibits a very complex motion, the ball can bounce between the two planes indefinitely without escape, and its motion is permanently chaotic [START_REF] Lehtihet | Numerical study of a billiard in a gravitational field[END_REF][START_REF] Tél | Chaotic dynamics: An introduction based on classical mechanics[END_REF]. However, in our case the complicated motion is observed for finite time, after that the ball "escapes" in either direction without further reversals. In this respect, our model is similar to a scattering process [START_REF] Bleher | Routes to chaotic scattering[END_REF][START_REF] Ott | Chaotic scattering: An introduction[END_REF]. Note, however, that despite the complexity emerging from the hither and thither motion, the A c c e p t e d m a n u s c r i p t behaviour does not become chaotic. In regard of this, the fly-wheel model is very similar to the infinite sequence of scattering disks investigated by Jung et al. [START_REF] Jung | Scattering one step from chaos[END_REF], where the existence of infinite periodic orbits does not imply chaotic scattering.

The presented fly-wheel model, although simplifies the motion of real balls, grasps the essential properties of bouncing balls. Our model could be enhanced, for example, by identifying the exact relation that determines the ratio ν of the outer to the inner moment of inertia. It is expected to depend on the elastic properties of the ball, but our model could be fine-tuned by introducing ν as a new, time-dependent variable. For example, one can speculate that this ratio might depend on the velocity and the spin of the impact, or it might even vary throughout one impact as a consequence of the varying normal force. Also, it is possible that the model could be enhanced by introducing more than two disks. However, the goal here was to find the simplest possible model that can perform the hither and thither motion.

Figure 2 :

 2 Figure2: Model of the fly-wheel. The two disks, separated in the figure by an annular gap for better visibility, can rotate with respect to each other around their common axis located in their centre, and they are attached to each other by torsional springs. The radius of the outer disk is R o , that of the inner is R i . The corresponding moments of inertia with respect to the axis are denoted by Θ o and Θ i . The positive direction of the rotations ϕ i and ϕ o is also indicated.

Figure 3 :

 3 Figure3: Free-body diagram during collision of the ball with the ground. The internal force A, maintaining the common motion of the disks, acts at the axis, and the internal moment M acts due to the torsional springs between the two disks. Two time-dependent forces act between the ground and the outer disk: the compressive force N and the friction force F .

Figure 4 :

 4 Figure 4:The deformation h as a function of time during one collision is shown with solid line, as described by Eq. (16). The time is measured in units of the total time τ of the collision, the deformation is measured relative to the maximum deformation h max . The approximation according to Eq. (17) is also shown with dotted line.

Figure 5 :

 5 Figure5: Normal force N (red line) and tangential force F (blue line) during an impact, as obtained from the fly-wheel model. The parameter values and the initial conditions are listed in the text. After an initial sliding phase (up to around 0.001 s) the tangential force is determined by the moment exerted by the internal rotation of the fly-wheel, which finally initiates sliding in the opposite direction after around 0.003 s. The discontinuity in the tangential force at around 0.001 s is the result of the jump between the sliding and rolling regimes.

Figure 6 :

 6 Figure6: The dimensionless tangential velocity β f of the contact point after the impact as a function of that before the impact (β c ). The parameters used for the simulation are the same as those used for Fig.5, but the velocity before the impact was always 4 m/s. For small incident angles (measured between the direction of the incoming velocity and the surface normal, β c 7), the ball starts to roll during the impact, for large angles (β c 7) it slides throughout the contact.

  us measure all distances in units of Θ o /R o m by introducing the dimensionless coordinates x = (R o m/Θ o )x, ỹ = (R o m/Θ o )y, and distances h = (R o m/Θ o )h. Also, we introduce dimensionless time by t = k/Θ o t, forces by F = F R o /k and Ñ = N R o /k, and parameters by ν

A c c e p t e d m a n u s c r i p t ( 21 )

 21 which gives the values () c of the variables right before the collision as a function of their values () f before the flight. The dimensionless duration of the impact is then computed from (15) as

  sliding occurs, and the dimensionless form of the equations of motion (5) become ẍ = -μN sgn(ζ ẋ + φo ), φ = -(1 + ν)φ -μN sgn(ζ ẋ + φo ), φo = -μN sgn(ζ ẋ + φo )φ, φi = νφ.

1 2 mR 2 0

 2 we would have ζ = 1 2 ν/(1 + ν).) This way, four independent parameters remain. In what follows, we set these parameters to the values μ = 0.8, α = 0.09246, ξ = 764.144, ν = 0.25,

Figure 7 :

 7 Figure 7: An example for the hither and thither motion of the fly-wheel. The fly-wheel changed the direction of its horizontal motion in 21 impacts. The parameter values were set according to (31), the initial conditions were those listed in the text, and v x0 = 1.0, ω 0 = 0.2.

Figure 8 :

 8 Figure 8: Dependence of the behaviour on the initial conditions. The colour coding, shown in b), indicates how many times the fly-wheel, released with the corresponding initial conditions ẋ(0) = v x0 and φo (0) = φi (0) = ω 0 , changes its horizontal velocity in the impacts before hopping off in one direction. The parameter values are as in[START_REF] Landau | Theory of elasticity[END_REF], further initial conditions are x 0 = 0, ẏ(0) = 0, y 0 = 543.48, and 500 × 500 initial conditions have been followed for a).

Figure 9 :

 9 Figure 9: Dependence of the behaviour on the initial spin ω 0 and on the parameter ν. The colour coding, shown in b), indicates how many times the fly-wheel, released with the corresponding initial condition φo (0) = φi (0) = ω 0 and parameter value ν rebounces before hopping off in one direction. The parameter values are those of (31) (except for ν and ζ), further initial conditions are x 0 = 0, ẋ(0) = 0, y 0 = 543.48, ẏ(0) = 0, and 500 × 500 initial conditions have been followed.

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t

Acknowledgements

Useful discussions with Tamás Tél and Róbert Németh are hereby gratefully acknowledged. We thank the financial support of OTKA under grant number K 068415.