Zhen Zhao

Caishan Liu

Bernard Brogliato

Planar dynamics of a rigid body system with frictional impacts. II. Qualitative analysis and numerical simulations

Keywords: multiple impacts, Coulomb's friction, experiments, numerical simulation

The objective of this paper is to implement and test the theory presented in a companion paper for the non-smooth dynamics exhibited in a bouncing dimer. Our approach revolves around the use of rigid body dynamics theory combined with constraint equations from the Coulomb's frictional law and the complementarity condition to identify the contact status of each contacting point. A set of impulsive differential equations based on Darboux-Keller shock dynamics is established that can deal with the complex behaviours involved in multiple collisions, such as the frictional effects, the local dissipation of energy at each contact point, and the dispersion of energy among various contact points. The paper will revisit the experimental phenomena found in Dorbolo et al. (Dorbolo et al. 2005 Phys. Rev. Lett. 95, 044101), and then present a qualitative analysis based on the theory proposed in part I. The value of the static coefficient of friction between the plate and the dimer is successfully estimated, and found to be responsible for the formation of the drift motion of the bouncing dimer. Plenty of numerical simulations are carried out, and precise agreements are obtained by the comparisons with the experimental results.

Introduction

The non-smooth dynamics exhibited in the dimer bouncing on a vibrated plate includes single and double impacts with friction, rolling contacts with or without slip, etc. In this paper, the numerical simulation for the bouncing dimer is implemented by using the theory developed in part I. Comparisons between the numerical and experimental results presented in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] will be carried out to illustrate and validate the scheme.

Basically speaking, the dynamics of the bouncing dimer will be divided into two parts: contact dynamics and impact dynamics. Both of them share a common set of dynamical equations that are constrained by different physical laws in both the normal and the tangential directions at each contact point. The Coulomb's friction law will be applied in both contact and impact dynamics to restrict the tangential motion at every contact point, while the normal constraints between them will differ owing to the dissipation and dispersion of energy in impact processes [START_REF] Ivanov | The problem of constrainted impact[END_REF][START_REF] Stronge | Impact mechanics[END_REF][START_REF] Ceanga | A new look at an old problem: Newton's cradle[END_REF][START_REF] Liu | Frictionless multiple impacts in multibody systems. I. Theoretical framework[END_REF][START_REF] Liu | Frictionless multiple impacts in multibody systems. II. Numerical algorithm and simulation results[END_REF]. For contact dynamics, the assumption of rigidity at contact surfaces permits a complementarity condition that has to be satisfied by the normal contact force and the normal acceleration at each contact point [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Pfeiffer | Multibody dynamics with unilateral contacts Wiley Series in Nonlinear Science[END_REF][START_REF] Brogliato | Nonsmooth mechanics[END_REF][START_REF] Glocker | Set-valued force laws[END_REF][START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF]. For impact dynamics, however, this complementarity condition will be replaced by the energetic coefficient of restitution to consider the local dissipation of energy, and a distributing law to consider the redistribution of energy among each contact point [START_REF] Liu | Frictionless multiple impacts in multibody systems. I. Theoretical framework[END_REF]. Additionally, the contact dynamics takes a set of ordinary differential equations to describe the global motion of rigid bodies in macro-size scales. The small time-scale at impacts allows one to introduce assumptions to ignore the effects from non-impulsive forces and the small variation of the configuration. Thus a set of differential impulsive equations can be established by using Darboux-Keller shock dynamics [START_REF] Keller | Impact with friction[END_REF][START_REF] Zhao | The analysis and simulation for three-dimensional impact with friction[END_REF].

The difficulty arising in Coulomb's friction is the indeterminacy of the tangential constraint when the relative tangential velocity vanishes [START_REF] Stronge | Swerve during three-dimensional impact of rough rigid bodies[END_REF][START_REF] Bhatt | Partitioning the parameter space according to different behaviors during three-dimensional impacts[END_REF][START_REF] Batlle | The sliding velocity flow of rough collisions in multibody systems[END_REF]. In order to identify the contact status at each contact point, local dynamical equations are established by transferring the governing equations into the local contact coordinate frame, in which the normal and tangential accelerations are separated and are directly connected with the contact forces. Once the relative tangential velocity disappears at a contact point, a correlative coefficient ðm st i Þ of friction for the contact point can be obtained by using the local dynamical equations with an assumption for a stick mode. If the correlative coefficient is located in the interior of the local frictional cone, a stick mode truly occurs when the tangential velocity vanishes. Otherwise, a reverse slip may appear, and the tangential constraint should be replaced by the Coulomb's friction law for slip. Once the tangential constraint is correctly identified, for the contact dynamics, the complementarity condition combined with the Coulomb's friction law will provide a solution for the local dynamical equations to obtain the normal and tangential contact forces that can make the numerical simulation advanced.

The process of impact is dominated by the contact point with maximal potential energy, which is linked with the 'time-like' independent variable [START_REF] Liu | Frictionless multiple impacts in multibody systems. I. Theoretical framework[END_REF][START_REF] Zhao | Experimental investigation of the Painleve ´paradox in a robotic system[END_REF]. Similar to the contact dynamics, the tangential constraint at the instant of the relative tangential velocity vanishing can also be identified by the correlative coefficient ðm st i Þ obtained from the local dynamical equations. If m st i is located in the interior of the frictional cone, the stick mode will appear in the contact point. Otherwise, the slip may reverse its direction when the relative tangential velocity disappears. Once the tangential constraints are determined, the combination of the distributing relationship and the energetic coefficient of restitution as well as Coulomb's friction law can make the impulsive differential equations be linked with a unique time-like independent normal impulse. Thus, the impulsive differential equations can be solved and output the post-impact velocities.

Owing to the special configuration of the dimer, a singularity will appear when it lies on the vibrated plate or experiences double impacts with two contact points sticking on the plate. This gives rise to a superstatic problem where there exists no solution in rigid body dynamics [START_REF] Brogliato | Numerical simulation of finite dimensional multibody nonsmooth mechanical systems[END_REF]. The other problem in simulation is how to identify a transition from a sequence of single impacts into a contact phase, which is associated with a so-called collisional singularity [START_REF] Falcon | Behavior of one inelastic ball bouncing repeatedly of the ground[END_REF]. These two problems are also discussed in detail in part I, and a scheme for dealing with them within the framework of rigid body dynamics is developed.

The paper is organized as follows. In §2, we present a brief summary for the basic equations of the bouncing dimer. In terms of the dynamical equations of the bouncing dimer and the correlative coefficient defined for stick mode as well as the experimental findings in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF], the static coefficient of friction between the dimer and the oscillated plate is estimated in §3. In §4, we present a qualitative analysis for the origins of the formation of the positive and negative drift motions in the bouncing dimer, in which the double impacts and the stick mode at contact points are found to play a significant role for the complex behaviour of the dimer. The comparison between numerical and experimental results is carried out in §5, and the difference for the slip coefficient of friction between impacts and contacts is claimed. Moreover, the influence of the initial and driving conditions on the formation of the persistent motions in the dimer is also investigated in §6, and some chaotic behaviours can be observed. Summaries and conclusions are given in §7.

The basic equations of the bouncing dimer

The dimer is shown in figure 1, in which two steel spheres are rigidly connected by a light glass rod and placed on a plate that vibrates with a sinusoidal waveform y p ðtÞZ A cosðutC aÞ. A and f are the amplitude and frequency of the plate oscillation, respectively. a is the initial phase angle of the plate oscillation related to the reference time for simulation. Let us denote (m b ,r), the mass of the sphere and its radius, respectively. The mass of the rod is m r and its length is (lK2r), which is varied to change the aspect ratio A r Z ð1C 0:5l=rÞ [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF]. J b Z 2m b r 2 =5, J r Z m r ðl K2rÞ 2 =12 are the inertias of the ball and the rod with respect to their mass centres, respectively. According to the theory developed in part I, the governing equations of the bouncing dimer are

M ðqÞ € q Khðq; _ qÞ Z K T ðqÞF C Q g ðqÞ; ð2:1Þ
where M(q) is a mass matrix; q Z ½x 1 ; y 1 ; q T ; (x 1 ,y 1 ) is the coordinate of the mass centres of the left ball; and q is the tilt angle of the dimer. K T (q) is the Jacobian matrix connecting the cartesian coordinates of the two potential contact points and generalized coordinates. hðq; _ qÞ and Q g (q) are the vectors of the inertial forces and the generalized forces due to gravity, respectively.

The local dynamical equations are established by separately expressing the normal and tangential accelerations at each potential contact point:

€ x n Z AF n C BF t C H n ð2:2Þ and € x t Z CF n C DF t C H t ; ð2:3Þ where € x n Z ½ € x n 1 ; € x n 2 T , € x t Z ½ € x t 1 ; € x t 2 T , F n Z ½F n 1 ; F n 2 T , F t Z ½F t 1 ; F t 2 T , KM K1 K T Z A B C D ! : H n Z ½H n 1 ; H n 2 T and H t Z ½H t 1 ; H t 2 T
are the terms of the normal and tangential accelerations induced by the applied and inertial forces.

For the contact dynamics, the local dynamical equations will be constrained by the following complementarity condition in the normal direction of every potential contact point:

€ d n $F n i Z 0; € d i R 0; F n i R 0; i Z 1; 2; ð2:4Þ
where d i Z x n i K y p ðtÞ are the gap functions for the relative displacements between the tips and the plate.

In the tangential directions of each contact point, the local dynamical equations will be constrained by the following Coulomb's law:

F t i ZKm i F n i sign v t i j v t i j ! ; if v t i s0; i Z 1; 2; if j F t i j % m s i F n i ; then v t i Z 0 and _ v t i Z 0; F t i ZKm i F n i sign _ v t i j _ v t i j ! ; if v t i Z 0 and _ v t i s0; i Z 1; 2; 9 > > > > > > > > > > > = > > > > > > > > > > > ; ð2:5Þ
where v t i is the relative tangential velocity, and v t i Z _ x t i , (iZ1, 2). m i O0 is the slip friction coefficient at the tip i. m s i O m i is the static coefficient of friction.

The local impact dynamics at the impulse level are expressed as

d _ x n Z A$dP n C B$dP t ð2:6Þ and d _ x t Z C$dP n C D$dP t ; ð2:7Þ w h e r e d _ x n Z ½d _ x n 1 ; d _ x n 2 T , d _ x t Z ½d _ x t 1 ; d _ x t 2 T , dP n Z ½dP n 1 ; dP n 2 T , dP t Z ½dP t 1 ; dP t 2 T , d _ x Z _ xðt iC1 ÞK _ xðt i Þ and dP Z ð t iC1 t i F dt Z ð t iC1 0 F dtK ð t i 0 F dt Z Pðt iC1 ÞKPðt i Þ ð 2:8Þ
are the changes of the relative velocities and normal impulses on [t i ,t iC1 ], respectively.

The local dissipation of energy at each contact point will be confined by the energetic coefficient of restitution,

e 2 s;i ZK W r;i E 0;i C W c;j ZK Ð P n i ðt f Þ P n i ðt c Þ _ d i dP n i E 0;i C Ð P n i ðt c Þ 0 _ d i dP n i ; ð2:9Þ
where W c,i %0 and W r,i R0 are the works done by the normal contact force at point j during the compression phase [0,t c ] and the expansion phase [t c ,t f ], respectively. E 0,i !0 is the initial potential energy at the contact point i.

The dispersion of energy during the impacts will be governed by a distributing law. For the two tips in the dimer with the same Hertz contact, the increment of the normal impulse, dP n 1 , at the left tip is connected to dP n 2 at the right tip by the following relationship: dP n i dP n j Z ðE ij ðP n i ; P n j ÞÞ 5=3 ; i Z 1; 2; i sj; ð2:10Þ

where E ij ZE i /E j is the ratio of the potential energies between contact points i and j. E i ðP n i ðtÞÞ is the residual potential energy at any instant t associated with a normal impulse P n i ðtÞ for the contact point i.

E i ðP n i ðtÞÞ Z E 0;i C 1 Tra ð P n i ðtÞ 0 _ d i dP n i ; ð2:11Þ
where Tra is a parameter to transfer the work done by the normal impulse into the potential energy: TraZ1 for a compressional phase ð _ d i ! 0Þ, and TraZ e 2 i for an expansion phase ð _ d i O 0Þ. The constraint in the tangential direction during an impact event is still governed by the Coulomb's friction law in equation (2.5), in which the quantities of contact forces are replaced by the corresponding increments of impulses.

Combining the distributing law with the impulsive differential equation and Coulomb's friction law makes frictional multiple impacts solvable. For the selection of the independent time-like variable, a guideline can be found in [START_REF] Liu | Frictionless multiple impacts in multibody systems. I. Theoretical framework[END_REF], in which the event of multiple impacts is dominated by a primary contact point that corresponds to the point taking the maximal potential energy among the various contacts.

Estimation of the static coefficient of friction

When a tip has a zero tangential velocity in both contact and impact phases, the tangential constraint depends on the property of the friction cone at the contact point (which is associated with the static coefficient of friction), and can be identified by using a correlation coefficient of friction. Let us first consider the correlative coefficient of friction in the case of the dimer with one contact point at i.

From equation (2.3) by setting F n j Z F t j Z 0 (isj ) and € x t i Z 0, we can express the correlative coefficient of friction m i;C st as

m i;C st Z j F t i j F n i Z 1 b 1 C b 4 C b 5 j b 6 KbH t i =F n i j ; ð3:1Þ where bZ m b l 2 C 4J b C 2J r , b 1 Z b=ð2m b C m r Þ, and b 4 Z ðm b C 0:5m r Þl 2 sin 2 q 2m b C m r ; b 5 Z 2rl sin q C 2r 2 ;
b 6 Z ð0:5l 2 sin q C rlÞcos q:

9 > = > ; ð3:2Þ
Since the value of F n i can be calculated by using the normal constraint of the complementarity condition, we can obtain the value of m i;C st and then compare it with the value of the static coefficient of friction m s i to identify the tangential status at the contact point i. Similarly, we denote the correlative coefficient of friction in the case of contact point i with a single impact as m i;I st , which can be explicitly expressed as

m i;I st Z j dP t i j dP n i Z b 6 b 1 C b 4 C b 5 : ð3:3Þ
In the case of double impacts in which the bouncing sphere collides against the oscillated plate while the other ball stays on the plate, the correlative coefficient of friction in this case is denoted as m DI st , which is

m DI st Z j dP t 1 C dP t 2 j dP n 1 C dP n 2 Z rl b 1 C 2r 2 j 1KðE 12 Þ 5=3 j ðE 12 Þ 5=3 C 1 : ð3:4Þ
When the dimer moves with a drift mode, the term H t i is related to _ q 2 that is very small and thus has little effect on the value of m i;C st . During the double impacts, the value of E 12 is also very small since one ball lies on the oscillated plate when double impacts occur. So we can approximately adopt a unified value of m st Z m i;I st to identify the property of the tangential constraint, which is only associated with the configuration of the dimer (the tilt angle q) and the geometrical parameter A r . By setting the steel density r steel Z7.8!10 3 kg m K3 and the glass density r glass Z1.0!10 3 kg m K3 , figure 2 shows the evolution of m st with A r and the angle q.

Now, let us estimate the stick coefficient of friction according to the experimental findings in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] and the unified correlative coefficient of friction m st . The experiments have shown that the dimer is directed towards the bouncing end if A r !5, otherwise the dimer moves in the backward direction. As shown in the experimental curves (see Dorbolo et al. 2005, fig. 3b), the reason for this is due to the different behaviour of the double impacts for small and large A r . In the case A r Z3.9, the horizonal velocity at the tip in the left ball can change with a positive value. This implies that a reverse slip occurs when the tangential velocity vanishes. However, for the case A r Z5.7 the positive value of the tangential velocity is null, thus a stick mode exists in the double impacts for the dimers with large A r .

If A r Z5 is set as the critical value for the occurrence of the stick mode, it is obvious that the drift behaviour of the dimer should also be changed before and after this value. Therefore, it seems reasonable that the value of m st at A r Z5 and qZ0 should be assigned on the stick coefficient of friction m s since the drift direction is changed after A r O5. From equation (3.3) and the knowledge that the slip coefficient of friction is near to the half value of m s , we estimate that the coefficients of friction approximately take the following values:

m s zm st ðA r Z 5; q Z 0Þ z0:24; m Z m st =2 Z 0:12: ð3: 5Þ
4. Qualitative analysis of the drift mode Now, let us qualitatively analyse the experimental phenomena shown in fig. 3b in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] for the tangential velocities of the dimer with different aspect ratios. In the case of A r Z3.9, the slip constraint of Coulomb's friction will be kept in double impacts. Therefore, the horizonal velocity can suddenly vary from a negative value to a positive one at double impacts. After that, the sequence of single impacts with the initial tangential velocities _ x K t 1 O 0 will appear in the motion of the dimer. These single impacts will be constrained by the tangential constraint, dP t 1 ZKm dP n 1 , and thus the tangential impulsive equations can be explicitly written as

d _ x t 1 ZK 1 b ½ðb 1 C b 4 C b 5 Þm C b 6 dP n 1 : ð4:1Þ
Since the coefficient before dP n 1 is always negative, the tangential velocity _ x t 1 decreases during the impact initiated from _

x K t 1 O 0. If the event _

x t 1 Z 0 is checked during the process of the impact, the slip constraint of Coulomb's friction will be kept, and a negative slip will appear in the contact point.

For the case of the dimer with A r O5, the stick constraint of Coulomb's friction will be kept in double impacts. Therefore, the value of the horizontal velocity is null after the event of the double impacts. After that, the sequence of single impacts with the initial tangential velocities _

x K t 1 O 0 will also appear in the motion of the dimer, and the tangential velocity _

x t 1 will decrease during the impact initiated from _

x K t 1 O 0. If the event _

x t 1 Z 0 occurs when the tilt angle q is small, the correlative coefficient will be smaller than the stick coefficient of friction, and thus the tangential constraint will be transferred into a stick mode that will make the tangential velocity at the contact point vanish when the impact finishes. Therefore, the horizonal velocity in the case of large A r seems to be null for a relatively long time. If the event _

x t 1 Z 0 appears in the dimer with a large tilt angle, the correlative coefficient may be greater than the stick coefficient, thus a slip mode will appear in the dimer with A r O5, and then a negative horizonal speed is allowed after a single impact.

If the single impact occurs when _ x K t 1 ! 0, the tangential constraint at the contact point becomes dP t 1 Z m dP n 1 , which will lead to the following tangential impulsive differential equation:

d _ x t 1 ZK 1 b ½ðb 1 C b 4 C b 5 ÞmK b 6 dP n 1 : ð4:2Þ
Obviously, the sign of the coefficient before dP n 1 will be influenced by the configuration of the dimer and the slip coefficient of friction.

If ðb 1 C b 4 C b 5 Þm! b 6 , the tangential acceleration of d _
x t 1 is guaranteed to be greater than zero, thus may produce a discrete event _

x t 1 Z 0 for the impact initiated from _ x K t 1 ! 0 and the tangential constraint will be changed during impacts. However, if ðb 1 C b 4 C b 5 ÞmO b 6 , the tangential acceleration of _

x t 1 will be always less than zero during the process of the impact, thus no variation of the tangential constraint occurs for the impact initiated from _

x K t 1 ! 0. In order to qualitatively analyse the property of the tangential motion in the impact initiated from _

x K t 1 ! 0, we can simplify the coefficient before dP n 1 by considering the dimer with a small tilt angle during the drift motion and ignoring the effects of the glass rod. Thus, we can set sin qz0, cos qz1 and m r Z0, and equation (4.2) becomes

d _ x t 1 Z f ðB r Þ 2ðB 2 r C 1:6Þ dP n 1 ; ð4:3Þ
where B r Zl/r, and the function f(B r ),

f ðB r Þ Z mB 2 r K2B r C 5:6m: ð4:4Þ
The sign of the function f($) depends on the roots of the equation f(B r )Z0.

Obviously, if mO 1= ffiffiffiffiffiffi ffi 5:6 p , no real root exists, so d _ x t 1 O 0 is guaranteed that will make _

x t 1 / 0 during the single impact initiated from _

x K t 1 ! 0. If 0! m! 1=
ffiffiffiffiffiffi ffi 5:6 p , there are two real roots:

B r 1 Z 2K ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1K5:6m 2 p 2m and B r 2 Z 2 C ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1K5:6m 2 p 2m : ð4:5Þ When B r 2 ½B r 1 ; B r 2 , d _
x t 1 ! 0 and thus the amplitude of _ x t 1 will increase during the impact with an initial condition of _

x K t 1 ! 0. So the discrete event of _ x t 1 Z 0 cannot appear in the process of the single impact. If B r ! B r 1 or B r O B r 2 , the tangential velocity _

x t 1 will approach zero, thus a transition of the tangential constraint from a slip mode into a stick mode is possible during the single impact. Whether the stick mode occurs or not depends on the static coefficient of friction m s , the value of A r and the tilt angle q.

From the estimation of the slip coefficient of friction, mZ 0:12! 1= ffiffiffiffiffiffi ffi 5:6 p and B r 2 ½B r 1 ; B r 2 . This illustrates that the single impact initiated from _

x K t 1 ! 0 cannot allow an event of _

x K t 1 Z 0, and the tangential constraint is not changed. Therefore, for the single impact initiated from _

x K t 1 ! 0, the amplitude of the tangential velocity in the tip of the dimer with small or large A r will increase until the impact finishes, since the event of _

x t 1 Z 0 cannot occur. A contact phase is confirmed by the detailed examination by [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] by measuring the electric resistance between the spheres and the plate, which is established by a sequence of single impacts. Obviously, the transition will occur if the intensity of the single impact is very small. During a (nonimpacting) contact phase, the tangential constraint at the contact point may also be changed due to the variation of the tangential velocity. For the dimer with a drift motion, the effects of inertial forces are in general very small, such that the qualitative property of the tangential constraint in contact phases is similar to the case of the dimer experiencing a single impact. Therefore, for a dimer with a small A r , the tangential velocity at the contact point will always decrease and can pass through zero to generate a reverse slip. However, for a dimer with a large A r , when the tilt angle is very small, the tangential velocity will be restricted to _

x t 1 R 0. Only when the tilt angle is large enough, the lock at the contact point can be released to permit a reverse slip.

Based on the above analysis, the behaviour of the tangential velocities shown in fig. 3b in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF], can be better understood. In the case A r Z3.9, no stick mode can appear in the double impacts or the sequential single impacts and the contact phases. So the tangential velocity can freely change without any restriction on its slip direction. Therefore, a positive tangential velocity can be obtained after double impacts, and will gradually decrease before its value equals zero, then increase its amplitude along a negative direction. Through a sequence of single impacts, a contact phase between the dimer and the plate can be established before or after the instant of _

x t 1 Z 0, and then a further cycle with a start of double impacts will begin: thus a positive drift motion is generated.

In the case A r Z5.7, however, the double impacts will make the tangential post-impact velocity vanish due to the occurrence of the stick mode. After that a sequence of single impacts occurs and then a contact phase can be established when the intensity of impacts becomes very small. In the situation of the dimer with a small tilt angle, the friction cone at the contact point only permits positive horizonal velocity, thus its value seems to be small or null for a relatively long time. Once the tilt angle is large enough, the restriction from the friction cone will be shaken off to admit a reverse slip. Then a negative slip can appear in the dimer.

Comparison between numerical and experimental results

According to [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF], the oscillated plate has a harmonic vibration with frequency f and amplitude A that is defined by a parameter GZ 4p 2 Af 2 =g, where g is the gravitational acceleration. In terms of the estimation of the coefficient of friction, and the restitution coefficient used in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] (eZ0.65 is the Newton's coefficient of restitution, and is identical to the energetic coefficient of restitution e s for a ball colliding against a fixed plate), we adopt the following parameters to carry out the numerical simulations: m s Z0.24, mZ0.12 and e s ZeZ0.65.

Let us first investigate the two cases exhibited in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] for the dimer with A r Z3.9 and 5.7. In both cases, the dimer is driven by the oscillated plate with the same parameters: fZ25 Hz and GZ0.9. The initial configuration of the dimer is denoted by a tilt angle q 0 and an initial height h 0 , the distance between the tip (the potential contact point) in the lower end of the dimer and the equilibrium position of the oscillated plate. The simulation starts at the instant when the waveform of the plate oscillation is expressed as y p (t)ZKA sin ut, where uZ1/f and A can be calculated according to the parameter G.

For the dimer with A r Z3.9, we compute l Z 2ðA r K1Þr Z 27:55 mm. Initially, the dimer is set with a configuration of q 0 Z58, and freely drop into the plate from a height h 0 Z0.25 mm. Figure 3 shows the numerical results for the vertical positions of the tips in the left (black solid curve) and right (black dashed curve) balls of the dimer. After several bounces to adjust the configuration of the dimer for participating impacts, a periodical vertical motion can be formed, in which the tip in the left ball of the dimer will stay on the plate, while the tip in the right ball will bounce off the plate periodically.

Except for the normal periodic motion, a novel horizontal drift can be observed in figure 4 for the tangential velocities at the left and right tips. From this figure, we can clearly observe the scenario for the complex motion of the dimer experiencing the drift mode: after the initial free motion of the dimer, the left sphere will collide repeatedly against the oscillated plate before the right ball reaches the plate. These single impacts will make the tangential velocities of both tips quickly decrease. After that, a strong collision between the right sphere and the plate occurs that makes the tangential velocities in both tips immediately change with a new positive direction. Then the right end of the dimer bounces far away from the plate, while the left end of the dimer experiences a sequence of single impacts. The next cycle for the right end of the dimer experiencing one collision per period of the plate oscillation begins.

Through the self-organization of the dimer motion experiencing a few cycles of the plate oscillation, the left ball will stay on the oscillated plate before the right ball collides against the plate. Then the single impact between the right ball and the plate will be transferred into double impacts. After that, a periodically complex motion will be formed that consists of double impacts, a sequence of single impacts at the left ball, and a contact phase between the left ball and the plate. Figure 3. The vertical positions of the left (solid curve) and right (dashed curve) contact points of the dimer with A r Z3.9 as a function of time obtained from numerical simulation, in which fZ25 Hz, GZ0.9, e s Z0.65, m s Z0.24 and mZ0.12. The thin broken curve is the vertical position of the plate.

The evolution of the tangential velocities at the tips in the left and right balls shown in figure 4 can also be better understood according to the qualitative analysis of the dimer dynamics. Since all the correlative coefficients in the dimer with A r Z3.9 are smaller than the stick coefficient of friction m s , the tangential velocities at the tips can freely pass through zero when the contact between the dimer and the plate is closed. So no stick mode appears in the tangential constraints of the contact points, such that the double impacts can make the tangential velocities of both tips positive, and a single impact at the left ball can make its horizonal velocity decrease rapidly, while the contact between the left sphere and the plate can slip reversely.

It is obvious that the stable persistent drift motion is generated by the repeated complex motion, which is periodic but rather subtle and much sophisticated in each cycle. Figure 5 shows the horizontal position of the left tip of the dimer in the case A r Z3.9. It is clear that the horizontal drift motion can be approximated as a uniform rectilinear motion with an approximately constant mean speed, huiz3.3 mm s K1 , obtained by dividing the distance x over 10 cycles by the corresponding time.

In comparison with the experimental results for the dimer with A r Z3.9 (see Dorbolo et al. 2005, fig. 3), a deviation at the quantitative level can be observed in the numerical results: the value of the mean drift speed obtained from simulation is lower than the one in experiments. From figure 4, one finds that the horizonal displacement is synthesized by the positive and negative slip of the left tip in each cycle of the complicated periodic motion. Therefore, the difference between them will determine the magnitude of the mean drift speed, which will be influenced by the amplitudes of the positive and the negative horizonal speed, and the durations for each part. Since no stick mode appears in the tangential constraints for the dimer with A r Z3.9, these quantities will be affected by the Figure 5. The horizontal position of the left tip of the dimer as a function of time in the drift mode for A r Z3.9 obtained from the simulation parameters: e s Z0.65, m s Z0.24, mZ0.12, fZ25 Hz and GZ0.9. slip coefficient of friction only. This illustrates that the estimated value mZ0.12 used in the simulation is not appropriate, and a modification for that will be discussed later.

Let us now investigate the case of the dimer with A r Z5.7, in which it will be directed to drift backwards from the bouncing end (a negative direction). The initial conditions and the parameters used in simulation are the same as for the case A r Z3.9: h 0 Z0.25 mm, q 0 Z58, lZ2(A r K1)rZ44.65 mm, m s Z0.24, mZ0.12 and e s Z0.65.

Similarly to the case A r Z3.9, a periodic bouncing motion in the normal direction is observed in figure 6, in which the bouncing end of the dimer will repeatedly collide against the oscillated plate, while the left sphere seems to stay on the plate.

However, the tangential motion of the dimer with A r Z5.7 will significantly differ from the one with A r Z3.9. Figure 7 depicts the evolution of the tangential velocities of the tips in the left and right balls for A r Z5.7. Let us check the complex motion of the dimer in this case.

The correlative coefficient m st with a configuration qZ58 is larger than the stick coefficient of friction. The initial single impact in the left sphere, after a free motion of the dimer, will therefore not take a stick mode in its tangential constraint. Consequently, the tangential velocity in the tip can change from zero to a negative value, and then the sequential detachment at the left ball will generate a sequence of single impacts that also makes the right ball approach the plate quickly. When the tilt angle q is so small to make the correlative coefficient m st be less than the stick coefficient of friction m s , a stick mode will appear in the single impacts at the tips of the left and right balls. Thus, the tangential post-impact velocities of the tips will vanish after the single impacts. Together with the decrease of the intensity in these single . The vertical positions of the left (solid curve) and right (dashed curve) contact points of the dimer with A r Z5.7 as a function of time obtained from numerical simulation, in which fZ25 Hz, and GZ0.9, e s Z0.65, m s Z0.24, mZ0.12. The thin broken curve is the vertical position of the plate. impacts, a contact between the left ball and the plate can be established. Then double impacts will appear in the dimer when the right ball reaches the plate again. Owing to the stick mode in the tangential constraint, the tangential velocities at both tips of the left and right balls will be equal to zero after the double impacts. The sequential single impacts at the left ball will make its tangential velocity approximately null due to the stick mode of the tangential constraint. Once the tilt angle q is large enough to make m st Om s , the stick mode at the tip of the left sphere will be shaken off, and the left ball can slip reversely. The repeated process will make the dimer drift backwards since the positive tangential velocity is absent in this case.

It is obvious that the complex periodic motions in the dimer with respective A r Z3.9 and 5.7 take a similar structure except for the variation of the tangential constraint during impacts and contact phases. Moreover, in comparison with the case A r Z3.9, the numerical simulations shown in figure 7 are close to the experimental results (see fig. 3 in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF].

Let us present an explanation for the deviation between the numerical simulations and the experimental findings for the dimer with the two different aspect ratios. For the dimer with A r Z5.7, the mean drift motion mainly depends on the negative amplitude of the tangential velocity at the tip of the left sphere before the double impacts occur. The negative tangential velocity is generated by the slip rolling motion that starts from a stick mode in the left sphere. The small error between the numerical and experimental results illustrates that the slip coefficient of friction for a contact phase used in simulation is approximately correct. However, for the dimer with A r Z3.9, the positive drift motion is synthesized by positive and negative slip motions that depend on the property of the slip friction in double impacts. Therefore, we may postulate the following for the frictional behaviour between bodies. Figure 7. The tangential velocities of the left (solid curve) and right (dashed curve) tips for the dimer with A r Z5.7 obtained from numerical simulation. The parameters for simulation are: e s Z0.65, m s Z0.24, mZ0.12, fZ25 Hz and GZ0.9.

Claim 5.1. Even though the stick coefficient of friction between bodies in impact and contact phases may take the same value, the slip coefficient of friction in impact phases should be less than the one in contact phases. 1In order to reflect the difference of the frictional behaviours in contact and impact phases, let us denote as m i the slip coefficient of friction during impacts, and m c the one in contacts. The slip coefficient of friction for contacts is kept as the value of m c Z0.12, and the one in impacts is estimated from a fitting method with a value of m i Z0.08, which can make the corresponding numerical simulation better coincide with the experimental results.

Under the same parameters as used in the above simulations except for the slip coefficient of friction during the impacts, figure 8 shows the horizonal velocities at the tips of the left ball for the dimers with the aspect ratios of A r Z3.9 and 5.7. Detailed observation for the horizonal velocities is presented in figure 9, zoomed from figure 8. Clearly, the new numerical results precisely agree with the experimental results in fig. 3a in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF].

The mean horizonal speed for the dimer with A r Z3.9 and 5.7 can be calculated by the horizonal displacements at the tip of the left ball shown in figure 10. In the case A r Z3.9, the mean horizonal speed is huiZ8.7 mm s K1 and the horizonal distance at the left tip over 120 ms is approximately x t 1 Z 1:044 mm. For the case Figure 8. The horizonal velocities at the tip of the left ball in the dimer with A r Z3.9 (solid curve) and 5.7 (dashed curve). The simulation parameters are: fZ25 Hz, GZ0.9, e s Z0.65, m s Z0.24, m c Z0.12 and m i Z0.08.

of A r Z3.9, huiZK6.7 mm s K1 approximately, and x t 1 ZK0:804 mm over a timer interval 120 ms. Both of them are very close to the values shown in fig. 3a in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF]. Some numerical simulations are also performed in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF], in which an event-driven algorithm is established by using Newton's equations of the dimer with the conditions for the occurrence of the various kinds of collisions 0.94 0.96 0.98 1.00 1.02 1.04 1.06 1.08 1.10 1.12 1.14 and the forces acting on the dimer. In particular, the transition between the stick and slip modes of the tangential constraint at contact points is also taken into account in their numerical schemes, such that the main feature of the drift motion in the dimer can be captured qualitatively. However, there is a significant deviation between their theory and experimental findings (see fig. 5b,c in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] because an oversimplified method is used to identify the transition of the tangential constraints, and to obtain the outcomes of the single and double impacts.

Let us check the influence of the aspect ratio A r and the frequency f on the dynamics of the dimer in drift mode after the modification of the slip coefficient of friction. The simulations are carried out for the dimers with various A r between 2.1 and 8.5 at three different frequencies of the plate vibration under the same energy GZ0.9. Each point shown in figure 11 is the scaled mean drift speed huið fg K1 Þ10 2 , where hui is the mean horizonal speed obtained by dividing the drift distance over ten cycles with the corresponding time. When the plate oscillates with a low frequency f and the dimer with a small aspect ratio A r , the excitation of the drift mode is more sensitive to the initial conditions. In particular, the drift mode cannot be excited for the dimer with A r !2.1 bouncing on the plate with fZ2.5 Hz.

A persistent stable drift mode strictly depends on the periodic behaviours of the dimer bouncing on the plate, in which each cycle of the repeated bouncing motion consists of complex motions such as double impacts, a sequence of single impacts and the contacts with friction effects. In general, a stable repeated cycle can be formed only if the dimer is released with an initial condition defined in a special range. If the initial conditions are not appropriate, the drift mode cannot be generated. At this stage, finding the basin of attraction of stable periodic motions is still an open problem. From figure 11, we can also find that the positive drift speed will decrease with the increase of the aspect ratio when A r 2(3.1,5), and the value of the frequency has little influence on the scaled mean drift motion. Meanwhile, an inflexion exists at A r Z3.1, in which the scaled mean drift speed in the case fZ25 Hz will rapidly decrease when A r changes from 3.1 to 2.9, and then increases by further diminishing the value of A r . In the cases of fZ50 and 75 Hz, the scaled mean drift speed seems to be a constant during the scope of A r 2(2.5,3.1), and then rapidly increases when the value of A r is further reduced.

In the case A r Z5, a transition from forward to reversal drift motion will suddenly occur. Indeed, the symmetric structure of the tangential velocity with small A r (see the relevant figures for A r Z3.9) collapses and is spontaneously transferred into the negative drift mode. Obviously, the symmetry breaking is due to the variation of the frictional behaviour at the contact point. In the negative drift mode, the absolute value of the scaled mean drift speed will also decrease when A r increases, and be significantly influenced by the value of the frequency of the plate. In the case fZ25 Hz, the absolute value of the scaled mean drift speed will decrease rapidly, and vanish when A r O6.2. This situation is mostly similar to a single pendulum, in which one end of the dimer is fixed on the oscillated plate, and the other swings beyond the plate and collides against it once every cycle. Figure 12 shows the normal motion for the dimer with A r Z6.5 released with zero velocities from an initial condition of h 0 Z2.25 mm and q 0 Z128, in which a small swing appears in the right ball even though the mean drift speed equals zero.

In the cases of fZ50 and 75 Hz, the occurrence for the zero mean drift speed are A r Z8.2 and 7.6, respectively. Under the same parameter G, the external and inertial forces applied in the dimer will vary with the frequency of the plate oscillation; thus, the critical value for the occurrence of the single pendulum mode will be translated. Figure 12. The vertical positions of the left (solid line) and right (dashed line) tips for the dimer with A r Z6.5 obtained from numerical simulation used in the following parameters: fZ25 Hz and GZ 0.9, e s Z0.65, m s Z0.24, m c Z0.12 and m i Z0.08. The thin broken curve is the vertical position of the plate.

Other stable persistent motion shown in the dimer

For a fixed A r , the mode of the persistent motion triggered in the dimer will be significantly influenced by the initial and driving conditions. For instance, the experiments in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] clearly show that the dimer dynamics is sensitive to the initial conditions, and the transition from a drift mode to the higher energy jump and flutter modes could be observed when the parameter G related to the plate oscillation is varied. In the following, we carried out some numerical simulations and made an attempt to give an illustration for the transition between different modes.

Once a stable persistent motion is established, the total energy involved in the dimer should be approximately constant such that a quasi-periodic motion holds. Let us analyse the energy involved in the drift, jump and flutter modes for the dimer.

In the case of the dimer with a stable drift mode, one ball stays on the plate, while the other ball swings beyond the plate and hits the plate once every cycle of the plate vibration. Most importantly, the swing motion of the dimer on the plate is seemingly periodic and the maximal height of the hitting end at every cycle seems to be kept approximately constant. As a simple way for the estimation of the total mechanical energies involved in drift mode we simplify the dimer dynamics as a system of an elastic ball bouncing off a fixed ground, in which the time interval between two sequential bounces of the 'elastic ball' is just equal to the period of the plate oscillation, T. The total mechanical energy in this simplified model is

E d Z mgh max Z mg$ 1 2 gðT=2Þ 2 Z 1 8 mg 2 T 2 ; ð6:1Þ
which can be thought of as the approximately constant energy involved in the dimer with a drift mode.

For the case of a jump mode, both the left and right balls will collide with the plate once every cycle of the plate oscillation. In this mode, the bouncing heights of the two balls are nearly the same. Since the glass rod is very light, the dimer resembles the situation of two separated particles bouncing on an oscillated plate with the same period that equals the period of the plate oscillation. From this simple model, the total mechanical energy of the dimer involved in the jump mode can be approximated as

E j Z 2mgh max Z 2mg$ 1 2 gðT=2Þ 2 Z 1 4 mg 2 T 2 : ð6:2Þ
The scenario for the dimer with a flutter mode is that both the left and right balls hit the plate out of phase at every other period of oscillation. In this mode the time interval between two sequential bounces is approximately equal to twice the period of the plate oscillation. Therefore, when the mode is simplified as the situation of two separated particles bouncing on an oscillated plate with the same period, the total mechanical energy of the dimer involved in the flutter mode can be estimated by

E f Z 2mgh max Z 2mg$ 1 2 gðT Þ 2 Z mg 2 T 2 : ð6:3Þ
Therefore, the ratio of the energy for the drift, jump, and flutter modes is approximately E d : E j : E f Z1 : 2 : 8, which is consistent with the statement in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF] and the numerical observation presented in this paper.

Obviously, the formation of a persistent motion means that the energy taken by the dimer will converge to a stable energy level that is maintained by its external environment. The convergence of the energy depends on the mechanism of the exchange between energies in the dimer and the oscillated plate. This exchange is implemented through the collisions between the dimer and the plate, and the stick-slip modes. Therefore, the initial conditions of the dimer, the plate oscillation and the property of the frictional collisions will significantly influence the formation of the persistent motion.

Let us first investigate the effects from the initial energy on the modes of the persistent motions in the dimer. For a dimer with A r Z3.9 under a plate oscillation with fZ50 Hz and GZ0.9, the numerical simulation shows that a stable drift mode can be formed when it is released with a small initial potential energy. However, the increase of the initial energy will trigger a transition from a drift to a jump mode. The vertical position for the dimer released with the initial conditions h 0 Z10.25 mm and qZ88 is shown in figure 13, in which a jump mode where both tips collide against the plate once every cycle appears in the dimer, and a horizonal reciprocating movement appears in the tips as shown in figure 14.

As noted in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF], the external energy involved in the plate oscillation also plays a significant role for the formation of the persistent motion. For example, Dorbolo's experiments have shown that, for a dimer with fixed A r Z3.9 and fZ25 Hz, a stable persistent drift motion can be observed over Figure 13. The vertical position of the tips in left (solid curve) and right (dashed curve) ball for the dimer (A r Z3.9, fZ50 Hz and GZ0.9) in a jump mode generated by an initial condition of h 0 Z10.25 mm and q 0 Z88. Parameters of simulations are e s Z0.65, m s Z0.24, m c Z0.12 and m i Z0.08. The thin broken curve is the vertical position of the plate. a range of G2[0.5,1.1], while the drift mode will spontaneously collapse to rest if G!0.5, or will be unstable with respect to the transition to the higher energy jump and flutter modes and their combinations if GO1.1. Figure 15. The horizonal position of the tip in the left (solid curve) and right (dashed curve) balls for dimer (A r Z3.9, fZ25 Hz and GZ0.2) generated by an initial condition of h 0 Z0.25 mm and q 0 Z128. Parameters of simulations are e s Z0.65, m s m c Z0.12 and m i Z0.08. The thin broken curve is the vertical position of the plate. Figure 14. The horizonal position of the tip in the left (solid curve) ball for dimer (A r Z3.9, fZ50 Hz, GZ0.9) in a jump mode generated by an initial condition of h 0 Z10.25 mm and q 0 Z88. Parameters of simulations are e s Z0.65, m s Z0.24, m c Z0.12, m i Z0.08.

Figure 15 shows the numerical results for a dimer with A r Z3.9 falling on the plate with an oscillation fZ25 Hz and GZ0.2. We can find that the dimer will quickly rest on the plate for the case of the plate with a weakening vibration.

Let us fix the values of A r Z3.9 and fZ25 Hz, but increase the magnitude of the plate vibration (GZ1.4). Figure 16 shows the simulation results for the dimer freely dropped on the plate with initial conditions h 0 Z0.25 mm and q 0 Z108. In this case, the vertical motion of the left ball seems to be harmonic with the plate oscillation, while the right ball will alternately collide against the plate with a big and a small bounce in each two adjacent cycles of the plate oscillation.

The above mode is unstable and sensitive to the initial conditions of the dimer. Let us carry out the numerical simulation by setting h 0 Z3.25 mm and q 0 Z108 for the same dimer dropped on the plate with the same oscillation. Figure 17 shows that a combination of jump and flutter modes appears in the dimer motions, in which the left ball collides against the plate once during every cycle, while the right ball bounces on the plate once during every two cycles of the plate oscillation. The increases of the initial energies involved in both the left and right balls will make the two adjacent bounces at the right ball be merged, and generate a big separation between the left ball and the oscillated plate.

The formation of an ordered persistent motion is also influenced by the initial value of the tilt angle in the dimer. A new mode is triggered when the dimer is released with a small height and a small tilt angle. Figure 18 shows the results obtained by setting the initial conditions with h 0 Z0.25 mm and q 0 Z58. The small initial tilt angle degrades the impact intensity at the right ball and thus makes the dimer's motion converge to a new stable status. The local picture, Figure 16. The horizonal position of the tip in the left (solid curve) and right (dashed curve) balls for dimer (A r Z3.9, fZ25 Hz and GZ1.4) generated by an initial condition of h 0 Z0.25 mm and q 0 Z108. Parameters of simulations are e s Z0.65, m s Z0.24, m c Z0.12 and m i Z0.08. The thin broken curve is the vertical position of the plate. figure 19 (zoomed from figure 18), shows that the dimer does not rest permanently on the plate, but leaves it with a small height during half of every cycle of the plate oscillation. Figure 17. The horizonal position of the tip in the left (solid curve) and right (dashed curve) balls for dimer (A r Z3.9, fZ25 Hz and GZ1.4) generated by an initial condition of h 0 Z3.25 mm and q 0 Z108. Parameters of simulations are e s Z0.65, m s Z0.24, m c Z0.12 and m i Z0.08. The thin broken curve is the vertical position of the plate. Figure 18. The horizonal position of the tip in the left (solid curve) and right (dashed curve) balls for dimer (A r Z3.9, fZ25 Hz and GZ1.4) generated by an initial condition of h 0 Z0.25 mm and q 0 Z58. Parameters of simulations are e s Z0.65, m s Z0.24, m c Z0.12 and m i Z0.08. The thin broken curve is the vertical position of the plate.

It is clear that the exchange between the energies involved in the dimer and the plate depends on the property of the impact between the dimer and the plate. Let us change the value of the coefficient of restitution e s . The dimer is freely dropped from a height h 0 Z0.25 mm and an angle q 0 Z48, with an oscillation fZ50 Hz and GZ0.9. The numerical simulation is carried out by setting e s Z0.8 Figure 20. The vertical position of the tips in left (solid curve) and right (dashed curve) of the dimer (A r Z3.9, fZ50 Hz and GZ0.9) in a flutter mode generated by an initial condition of hZ4.9 mm and q 0 Z48. The simulation parameters are e s Z0.8, m s Z0.24, m c Z0.12 and m i Z0.08. The thin broken curve is the vertical position of the plate. and keeping other simulation parameters with the same values as used in previous cases. Figure 20 shows the vertical positions of the tips in the dimer. Clearly, a flutter mode is generated, in which both the left and right balls hit the plate out of phase at every other period of oscillation.

Summary and conclusions

In terms of the theory presented in part I, the non-smooth dynamics of a dimer bouncing on a vibrated plate is qualitatively analysed and simulated in detail. The frictional property between the plate and the dimer is successfully estimated based on the experimental findings combined with the theoretical developments for the non-smooth dynamics, which leads to a precise agreement between numerical and experimental results. Theoretical analysis and numerical investigations further reinforce the statement presented in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF], in which the friction and multiple shocks play a prominent role for the persistent motions of the bouncing dimer. Moreover, chaotic behaviours are observed in numerical simulations, where periodical motions can be triggered in the dynamics of the dimer and are extremely sensitive to the initial and driving conditions, as well as the physical parameters related to the contact/ impact behaviour.

In comparison with the numerical scheme provided in [START_REF] Dorbolo | Dynamics of a bouncing dimer[END_REF], this paper provides a systematical theory that can correctly deal with the problem of the frictional multiple impacts. The precise agreements between numerical and experimental results indicate that the distributing law found in [START_REF] Liu | Frictionless multiple impacts in multibody systems. I. Theoretical framework[END_REF] for frictionless multiple impacts can be extended into the case with friction, and supports the multivalued part of Coulomb's graph is the intrinsic property of dry friction and should be carefully considered in simulations. Further research work will focus on the aspect of establishing a uniform framework for non-smooth dynamics, in which the many singularities induced by friction should be investigated carefully.

Figure 1 .

 1 Figure 1. A dimer bouncing on a vibrated plate.

Figure 2 .

 2 Figure 2. The value of m st with the variation of A r and the angle q.

Figure 4 .

 4 Figure4. The tangential velocities of the left (solid curve) and right (dashed curve) tips for the dimer with A r Z3.9 obtained from numerical simulation. The parameters for simulation are: e s Z0.65, m s Z0.24, mZ0.12, fZ25 Hz and GZ0.9.

  Figure6. The vertical positions of the left (solid curve) and right (dashed curve) contact points of the dimer with A r Z5.7 as a function of time obtained from numerical simulation, in which fZ25 Hz, and GZ0.9, e s Z0.65, m s Z0.24, mZ0.12. The thin broken curve is the vertical position of the plate.

Figure 9

 9 Figure 9. A local figure zoomed from figure8for the horizonal velocity at the tip of the left ball in the dimer with A r Z3.9 (solid curve) and 5.7 (dashed curve), respectively.

Figure 11 .

 11 Figure11. The scaled mean horizonal speed huið fg K1 Þ10 2 versus A r for three values of driving frequency. Parameters used in simulations are: GZ0.9, e s Z0.65, m s Z0.24, m c Z0.12 and m i Z0.08.

Figure 19

 19 Figure19. A local picture zoomed from figure18for the horizonal position of the tip in the left (solid curve) and right (the dashed curve) balls for dimer (A r Z3.9, fZ25 Hz and GZ1.4). The thin broken curve is the vertical position of the plate.

It is not easy to give rigorous support for this statement unless a lot of experiments are carried out. Intuitively, the strong interaction during impacts will weaken the tangential resistance due to the local plastic deformation, such that the slip coefficient of friction during impacts will be lessened.

We would like to express our appreciation to the anonymous reviewers for their beneficial comments and suggestions. The good experimental report and the hints from their numerical simulations by Dorbolo et al. in 2005 establish a solid foundation for the theoretical developments presented in this paper. The support of the National Science Foundation of China (10772002) and PHR (IHLB) is gratefully acknowledged.


Rotation Curves and Dark Matter

Keywords: Dark Matter, Rotation Curves, Gravitational Potential, Mass of galaxy. Pacs numbers: 04.60.-m, 04.65.+e, 11.15.-q, 11.30.Ly May 2017

In present paper we argue that to explain the shape of the Rotation Curves (RC) of galaxies, there is no need to involve the concept of dark matter. Rotation curves are completely determined by the distribution of baryon matter and gas kinetics. Such parameters of the galaxy as barion mass and its distribution can be easily calculated from the observed RC. We show the extended parts of RCs to be just a wind tails, formed by gas of the outer disks in assumption that it obeys the laws of gas kinetics. As examples, the Galaxy, NGC7331 and NGC3198 are considered. We calculate total mass of the Galaxy and find it to be 23.7x10(10)M sun. For the NGC7331 and NGC3198 the calculated total masses are 37.6x10(10)M sun and 7.7x10(10)M sun respectively. Consequences for cosmology are discussed.

Introduction

From the very beginning of its discovery in 1932 by Jaan Oort and confirmation made by Fritz Zwicky in 1933, the phenomenon of Dark Matter (DM) was widely discussed. Unfortunately, in the past 85 years, there has been no significant progress in understanding the nature of DM and the dynamics of galaxies. In present paper we have tried to fill this gap and offer a physically reasonable explanation for the DM phenomenon.

To begin with, we would like to mention here some arguments against existence of DM, at least in it's actual meaning.

-It is well known that the halo stars (the old population of the galaxy) move more slowly than the young population of the disk does. This fact clearly indicates an unsatisfactory explanation of the Rotation Curves (RC), as provoked by the spherically distributed DM.

-An immense dispersion of the star velocities beginning at distances of the order of 10 kpc, which can not be explained with a homogeneous, spherically distributed DM. Such a spread can be explained probably by the presence of a compact structure such as the arms of galaxy, which are baryonic by nature.

-Recently it was shown that there is a significant correlation between the features in the RC and the spiral structure of the baryonic disk. As the authors say: "The dark and baryonic mass are strongly coupled" [START_REF] Mcgaugh | The radial Acceleration Relation in rotationally supported galaxies[END_REF], [START_REF] Sancisi | Dark Matter in Galaxies[END_REF], [START_REF] Mcgaugh | [END_REF], [START_REF] Möller | The kinematical structure of gravitationally lensed arcs[END_REF] but it is impossible in the framework of the usual DM paradigm.

-This year the rotation curves for high redshifted galaxies were reported and it was clearly shown that a large fraction of massive high -redshifted galaxies are strongly baryon-dominated (see [START_REF] Genzel | Strongly baryon-dominated disk galaxies at the peak of galaxy formation ten billion years ago[END_REF] and references therein). These data clearly indicate the absence of the DM in the early epochs and its mystical appearance in our time.

-The FERMI experiment for search of anniquilation of DM and anti-DM, clearly show negative result [START_REF] Fermi | Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with Fermi-LAT[END_REF] All of the above clearly indicate the need to revise the concept of DM and a more accurate simulation of the galactic baryon component.

In present paper we prove that the RC of spiral galaxies are just the wind tails of the baryon gas following the preceding baryon matter, and we do not need the DM concept to explain rotation curves. We base our proof not only on numerical estimates, but also on exact calculations for galaxies and a comparison of the results with observations. The article is organized as follows:

In part 2 we discuss RC produced by spherically distributed matter, we suggest an analytical solution for the density function under condition V ⊥ = const, and compare it with those obtained with numerical simulations.

Part 3 is devoted to the construction of a realistic and convenient density function of the baryon component in a cylindrical coordinate system. We compare it with the Miyamoto-Nagai function by example of the Galaxy. Analytical solution for the RC, that corresponds to the density function is suggested. Expression for the total mass of S-type galaxies is obtained

The main part 4 is devoted to the gas dynamics and relationship of the baryon gas density distribution function to the RC. Obtained results are in excellent agreement with the observational data.

In conclusion the main results of this paper are summarized and important consequences are discussed.

Rotational curves in the case of spherical sym-metry

The case, when the density distribution has the spherical symmetry deserves a separate consideration because on the one hand the DM is thought to be spherically distributed, as it is supposed in all numerical simulations, and on the other hand we need it to compare with the density function obtained in part 4 of our paper. So we consider here the spherically symmetric case in details.

It is well known that for point -like source the gravitational potential can be written as

du = Gdm(r) r . (1)
Let's consider spherically symmetrical density σ(ρ) in [g • cm -3 ]. In this case the mass differential can be expressed as

dm = σ(ρ)ρ 2 sin θdρdθdφ . ( 2
)
For this reason we have the following potential function for a spherically symmetric system:

du = Gσ(ρ)ρ 2 sin θdρdθdφ R 2 + ρ 2 -2Rρ cos θ , ( 3
)
where R is the distance from the observer to the origin of the sphere under consideration.

Performing elementary integration with respect to the angles θ and φ, we obtain:

du = 4πGσ(ρ) R ρ 2 dρ, (4)
or

u = 4πG R R 0 σ(ρ)ρ 2 dρ, (5)
which gives us in particular case of constant density the following relation for the square velocity:

V 2 ⊥ = R du dR = 4πGσR 2 , (6)
Now let's consider how should density depends on distance to maintain condition V ⊥ = const. in the case of spherically distributed matter.

The DM density function (which could produce approximately constant RC in a certain range of distances) is widely discussed and applied to model DM haloes, see for example [START_REF] Navarro | The structure of cold dark matter halos[END_REF], [START_REF] Merritt | Empirical models for dark matter halos. I. nonparametric construction of density profiles and comparison with parametric models[END_REF], [START_REF] Katz | Testing feedback-modified dark matter haloes with galaxy rotation curves: estimation of halo parameters and consistency with Lambda-CDM scaling relations[END_REF] and references therein. Unfortunately to date, such density functions are calculated only from numerical simulations. We will not offer here an exhaustive description of all such functions, but mention just two most commonly recognized. The first were suggested by Navarro, Frenk and White [START_REF] Navarro | The structure of cold dark matter halos[END_REF] and has the form

σ N F W (r) = σ 0 r rs (1 + r rs ) 2 . (7)
Another function was obtained recently by Di Cintio et al. [10], [11] and was compared with the first one in [START_REF] Katz | Testing feedback-modified dark matter haloes with galaxy rotation curves: estimation of halo parameters and consistency with Lambda-CDM scaling relations[END_REF]. As it was stressed in [START_REF] Katz | Testing feedback-modified dark matter haloes with galaxy rotation curves: estimation of halo parameters and consistency with Lambda-CDM scaling relations[END_REF], the function of Di Cintio et al. much better approximates the observed RC. This density distribution of DM can be represented as:

σ DC14 (r) = σ 0 ( r rs ) |γ| (1 + ( r rs ) α ) β-γ α . (8)
In the case of DM dominated galaxy, this distribution is characterized by the following values of the parameters γ ≈ -0.88, α ≈ 1.4 , β-γ α ≈ 2.6 . It should be emphasized again that these functions were obtained from numerical simulations on the basis of the requirement to describe an approximately constant velocity in a strictly defined range of distances, in order to approximate RCs of real galaxies. Let's calculate analytically the density function which corresponds to constant velocity. By taking into account that

V 2 ⊥ (R) = R du dR , (9)
and requiring that V ⊥ = const , we immediately obtain from the eq. ( 5) the following differential equation for the density σ(r):

σ ′ (R)R 3 + 2σ(R)R 2 = V 2 4πG . (10)
Or by normalizing the variable for convenience of comparison x = r/r s , we have:

σ ′ (x)x 3 + 2σ(x)x 2 = V 2 4πGr 2 s . (11)
Integrating (11), we find

σ(x) = V 2 4πGr 2 s ln( r rs ) ( r rs ) 2 . ( 12
)
This is exact expression for the density function, which comply condition V ⊥ = const over the entire range of distances. As can be seen, expressions [START_REF] Navarro | The structure of cold dark matter halos[END_REF], [START_REF] Merritt | Empirical models for dark matter halos. I. nonparametric construction of density profiles and comparison with parametric models[END_REF] and [START_REF] Miyamoto | Three-Dimentional Models for the Distribution of Mass in Galaxies[END_REF] have similar behavior for the argument r rs greater than one.

Rotation curve in the case of cylindric symmetry

It is well known that the baryon component of galaxies consists of stars and gas.

The stars can be considered as collisionless gas, while the real gas also obeys the laws of gas dynamics and should be described by Fick's laws. Due to the fundamental differences, these two components should be considered separately and in this part of paper we will consider only the stellar component. Physics of the gas component will be considered in the 4-th part of the paper. Now consider a thin disk of radius a , characterized by thickness 2z 0 = 2b. Let dm be a point -like mass inside of the disk and r is the distance from the mass to observer (which also is located inside of the disk at distance R from the center of galaxy). Than the distance observer -point-like mass can be written as:

r 2 = z 2 + R 2 + ρ 2 -2Rρ cos ϕ, (13)
where ρ , z and ϕ are cylindrical coordinates of the mass under consideration.

For the convenience, we split the potential formed by the mass in the point of observer, into longitudinal and tangential components u = u + u ⊥ .

It is clear that tangential component will not affect on the RC. For longitudinal component we have

u = u R -ρ cos ϕ r . ( 14
)
But for the point -like mass the potential is u = Gm/r, so

du = Gdm R -ρ cos ϕ r 2 , (15)
or in it's complete form

u = G R 0 b -b 2π 0 σ(ρ, z)(R -ρ cos ϕ)ρdρdϕdz z 2 + R 2 + ρ 2 -2Rρ cos ϕ , (16)
and integrating over ϕ we obtain:

u = 2πG R R 0 b 0 σ(ρ, z)ρdρdz. (17)
In order to follow further, we need the density distribution function be defined. To the best of our knowledge, the most used one to describe observed galaxies, is the Miyamoto -Nagai density function [START_REF] Miyamoto | Three-Dimentional Models for the Distribution of Mass in Galaxies[END_REF]. Unfortunately, this function is not convenient for analytical calculations and therefore we approximate it by a factorized function of the same degree that does allow simple integration:

σ(ρ, z) = 10 10 M ⊙ (γt 2 + 1) 3/2 α k (β k x 2 + 1) 3/2 k . ( 18
)
Here γ, α k , β k are parameters of approximation and dimensionless variables are x = ρ/a and t = z/b. The density function of Miyamoto -Nagai (M-N) is widely used in calculations, and therefore it will be useful to compare our model [START_REF] Begeman | HI rotation curves of spiral galaxies[END_REF] with the well known M-N density function for the Galaxy. For the density function [START_REF] Begeman | HI rotation curves of spiral galaxies[END_REF], the integration over z can be carried out easily:

I = b 0 1 (γt 2 + 1) 3/2 dz = bt γt 2 + 1 1 |= 0 b √ γ + 1 , (19)
and we obtain

u = ηb R √ γ + 1 R 0 k α k (β k x 2 + 1) 3/2 ρdρ, (20)
where we introduce the constant η = 2πG10 10 M ⊙ . Integration over ρ gives

u = ηa 2 b R √ γ + 1 k α k β k 1 - 1 (β k R 2 a 2 + 1) 1/2 . ( 21
)
This is gravitational potential produced by galaxy in the point R. The required RC can be found now both from [START_REF] Kroupa | The failures of the standard model of cosmology require a new paradigm[END_REF] and from (21) as

V 2 ⊥ = ηa 2 b R √ γ + 1 k α k β k 1 - 3 2 β k R 2 a 2 + 1 (β k R 2 a 2 + 1) 3/2 , ( 22
)
where all distances are measured in kpc. This expression was obtained for the density function [START_REF] Begeman | HI rotation curves of spiral galaxies[END_REF] characterized by the parameters γ, α k and β k , and it can be used immediately to calculate RC. Now we introduce for convenience new coefficients α * k and β * k in which we include the parameters of the model. Namely let α * k = bα k / √ γ + 1 , and

β * k = β k /a 2 .
In this case [START_REF] Lipovka | Physics on the Adiabatically Changed Finslerian Manifold and Cosmology[END_REF] became:

V 2 ⊥ = η R k α * k β * k 1 - 3 2 β * k R 2 + 1 (β * k R 2 + 1) 3/2 . ( 22a
)
Plotted with relation (22a) RC for our Galaxy is shown in Fig. 2, in comparison with the observed RC [START_REF] Clemens | Massachusets-Stony Brook Galactic Plane CO Survey: The Galactic Disk Rotation Curve[END_REF], [START_REF] Blitz | Catalog of CO Radial Velocities Toward Galactic HII Regions[END_REF], [START_REF] Burton | Carbon Monoxide in the Galaxy. III. The Overall Nature of its Distribution in the Equatorial Plane[END_REF] As can be seen, the obtained RC perfectly coincides with the observational data in the range R < R 0 ≈ 15 kpc , while one can see the discrepancy for the outer part of the disk (We define R 0 as a distance, at which the tangential velocity

V ⊥ = V K if R < R 0 and V ⊥ = (V K +V d ) > V K if R > R 0 where V K is the Keplerian velocity).
It can be explained by the fact that the RCs for distances R > R 0 are measured mainly with HI 21cm. line and formed by gas. For this reason this part of the RC can not be obtained from simple gravitational model, and more sophisticated physics based on the gas kinetics should be applied. We consider it in details in the part 4, but now to conclude this part we would like to suggest evaluation for the total mass of our Galaxy. For the density distribution model used for better fit of the observed RC presented in Fig. 2 (solid line) we obtain, integrating [START_REF] Begeman | HI rotation curves of spiral galaxies[END_REF]:

M G = 4π10 10 M ⊙ k α * k β * k 1 - 1 β * k a 2 + 1 . ( 23
)
By substituting the parameters found above and taking for the Galaxy a = 15 kps (it should be noted here that parameter a weakly affects the calculated mass), we obtain the baryon disk mass (actually this is the total mass because the mass of the gas is negligible) of our Galaxy M G = 23•10 10 M ⊙ which perfectly agrees with the Oort limit.

Rotation curves measured with HI line 4.1 Preliminary estimates

To begin with it should be stressed again that there are two very different components of galaxy population: stars and gas which are used to measure the RC of galaxy in optics and in radio respectively. The first component is drived only by gravitation potential, whereas to describe the second one we should take into account collisions and the diffusion equations should be involved into consideration. Rough estimate of the mean free path time for a hydrogen atom t f p = (N σV t ) -1 (here N is the density of the gas in sm -3 , σ is cross-section for elastic collision and V t is the mean thermal velocity of the atom) gives t f p ≈ 1.3 • 10 10 /N (sec) = 4.1 • 10 2 /N (yrs). For the typical HI density outside of the R 25 : N ≈ 10 -3 -10 -4 we obtain evaluation t f p ≈ 4.1 • (10 5 -10 6 ) yrs. Thus, one can see that for description of the neutral hydrogen, located in the outer R > R 25 R 0 part of galaxy's disk, the complete gas dynamic equations should be used for correct description of the RC. Now let's made another estimation to answer the question : "Will the gas be able to follow the underlying falling baryon matter to form the wind tails?". From observations we know that baryon matter of a S-type galaxy moves along a spiral. In this case we may model it as a piston (baryon matter in inner R < R 0 part of galaxy) that moves in the spiral tunnel with ideal walls, and is followed by the HI gas (here we will not consider the processes of star formation, that dilute the gas component, but we only note here that accounting for such processes will increase the effects we are discussing). The mean acceleration of the "piston" for typical galaxy can be evaluated as < w >= ∆V /∆t = (200 Km/s) / (10 9 yrs) = 10 -9 (cm/s 2 ) . By taking into account the evaluation of t f p made before, we can estimate the variation of the piston's velocity ∆V during the mean free path time of the hydrogen atom. We have ∆V =< w > t f p = 10 -9 • 2 • 10 10 /N = 20/N (cm/s). For typical density N = 10 -3 (cm -3 ) we obtain ∆V = 2 • 10 4 (cm/s) << V t ≈ 10 6 (cm/s). So one can conclude that gas will follow the "piston" if the gas density is enough high: N (cm -3 ) 2 • 10 -5 • (10 6 /V t ) . These were rather crude assessments, suggested here to show simplistically the physics of processes. To conclude this part we would like to stress that in consequence with these simple evaluations, the gas, driven by collisions, can easily follow the underlying baryon matter. Actually the gas under consideration forms the wind tail which is rigidly follows the underlying baryon matter that is driven mainly by the gravity at distance R 0 . This way the absence of RC of S-type galaxies in early universe can be explained easily. Rough estimate of distance over which the wind tail (or, the same RC) can spread is t • V t ≈ 10 10 yrs • 3 • 10 7 • 10 6 ≈ 3 • 10 23 (cm) = 100 kpc. This trivial evaluation clearly shows why the RC measured with HI line are seen in our epoch, but can not be observed in early universe, when t < 10 10 yrs., as it was recently reported [START_REF] Genzel | Strongly baryon-dominated disk galaxies at the peak of galaxy formation ten billion years ago[END_REF].

Gas density as a function of distance for constant RC

Now we consider detailed description of the process by using the diffusion differential equations. From observations we know that gas in outer parts of spiral galaxy follows a spiral which we write as

R = R 25 e k(ϕ) , ( 24
)
where the distance R 25 has usual definition. In this case for the total longitudinal V and tangential V ⊥ velocity one can write V = k ′ V ⊥ , where k ′ is ∂k/∂ϕ . But the longitudinal velocity is given by the Fick's first law of diffusion:

V = - D N ∂N ∂R , (25)
or

k ′ V ⊥ = - D N ∂N ∂R , ( 25a
)
where D is the diffusion coefficient, N is density and R is distance in the cylindric coordinate system. Formally speaking one can put into equation (25a) the HI gas density, measured for different spirals and obtain RC, but it is methodologically more correct and easier to solve the inverse problem. Namely, we are interested in the question: "which HI column density function corresponds to the case of constant rotation curve of baryon matter in absence of DM for an S-type galaxy?" By taking into account the fact that D = const for very depleted gas and assuming that V ⊥ = V K0 = const., (here V K0 is Keplerian velocity at the distance R 25 ) we obtain

D n 1 R 25 ∂n ∂r = k ′ V K0 , (26)
where we introduce n = N/N 25 , r = R/R 25 , and by definition

V K0 = M G R 25 . ( 27
)
Integrating ( 26) we obtain:

n = n 0 exp - V K0 R 25 k ′ D r , (28)
where

n 0 = exp V K0 R 25 k ′ D . ( 29
)
As one can see this density distribution distinct dramatically of that, obtained for the spherically distributed DM [START_REF] Miyamoto | Three-Dimentional Models for the Distribution of Mass in Galaxies[END_REF] written for a collisionless DM gas. Now we are ready to find the column density formed by distribution (28).

Column density

By definition the column density is

N L = 2 L 0 N dl, ( 30
)
where N is density function. By taking into account that l 2 = R 2 -ρ 2 , the eq. ( 30) can be rewritten as:

N L = 2N 25 R 25 rmax r=ρ/R25 nrdr r 2 -ρ 2 R 2 25 . ( 31
)
Expanding the expression a row:

1 1 - ρ rR25 2 ≈ 1 + 1 2 ρ rR 25 2 + 3 8 ρ rR 25 4 , (32)
we obtain

N L (ρ) ≈ 2N 0 R 25 rmax r=ρ/R25 n(r) 1 + 1 2 ρ rR 25 2 + 3 8
ρ rR 25 4 + ... dr, (33)
where n(r) is given by relation (28). So, the column density can be expressed as

N L ≈ N L1 + N L2 + N L3 , (34)
where

N L1 = 2N 25 R 25 rmax r=ρ/R25
e -κ(r-1) dr,

N L2 = N 25 R 25 ρ R 25 2 rmax r=ρ/R25 (35)
e -κ(r-1) dr r 2 , (36)

N L3 = 3 4 N 25 R 25 ρ R 25 4 rmax r=ρ/R25 e -κ(r-1) dr r 4 , (37)
where

κ = V K0 R 25 k ′ D . (38)
these expressions can be integrated and we obtain for the first term:

N L1 ( ρ R 25 ) = 2N 25 R 25 κ e -κ( ρ R 25 -1) 1 -e -κ ρ R 25 ∆V V , (39)
where ∆V /V is the relative bandwidth by velocity.

Estimating the integrals (36) and (37), it is easy to show that the second, third and higher terms of our expansion, will have the same exponential behavior (∝ exp(-κ(ρ/R 25 -1))), but in absolute value they will be much smaller than the first term (39). For this reason for our aim we can take expression (39) as a good approximation for the column density. The column density (39) is compared with observed values [START_REF] Bigiel | Extremely inefficient star formation in the outer discs of nearby galaxies[END_REF] for NGC7331 (fig. 3) and NGC3198 (fig. 4).

These galaxies were chosen because they are placed edge on in respect to the observer. Due to this circumstance, in this case, there is no need to take into account the angle of inclination of galaxy, that simplifies the model.

As can be seen at fig. 3 and fig. 4, obtained density distribution excellently approximates the observed column densities for very different galaxies, characterized by different slopes of the column density function. So we can conclude that RC are formed by wind tails of falling gas that obeys the diffusion equations, and we do not need dark matter to explain the rotation curves of galaxies. Now we estimate the masses of two galaxies mentioned above from their measured RC suggested in [17], [START_REF] Begeman | HI rotation curves of spiral galaxies[END_REF] and [START_REF] De Blok | High-resolution rotation curves and galaxy mass models from[END_REF] by using previously obtained relation [START_REF] Lipovka | Planck Constant as Adiabatic Invariant Characterized by Hubble's and Cosmological Constants[END_REF]. The coefficients α * k and β * k we immediately find from approximation of the RCs for these two galaxies. Figures 5 and6 demonstrate results of such approximation for NGC7331 and NGC3198 respectively.

[M sun pc -2 ] Galactocentric Radius [R 25 ] 1,2 1,4 1,6 1,8 2,0 0,01 0,1 1 10 HI
Figure 3: Measured (squares) [START_REF] Bigiel | Extremely inefficient star formation in the outer discs of nearby galaxies[END_REF] and calculated with (39) (solid line) column density of HI for NGC7331.

The thick straight line in the figures shows wind tails formed by gas. As can be seen, wind tails extend exactly to the distance where the column density function has the exponential form (39) (see also fig. 3 and fig. 4).

Obtained coefficients for NGC7331 are α * 1 = 0.2317 , β * 1 = 0.112 , α * 2 = 6.358 , β * 2 = 28.8, and for NGC3198 we find α * 1 = 0.2317 , β * 1 = 0.112 , α * 2 = 6.358 , β * 2 = 28.8. Now the masses of these galaxies can be obtained immediately with relation [START_REF] Lipovka | Planck Constant as Adiabatic Invariant Characterized by Hubble's and Cosmological Constants[END_REF]. For NGC7331 we have M 7331 = 37.6 • 10 10 M ⊙ and for NGC3198 the total mass is M 3198 = 7.7 • 10 10 M ⊙ .

Conclusions

As it is known, the need for the DM arises from the three observable problems (we will not mention here the Cosmology and problem of the observed structure formation for the reasons stated below).

1,0 1,2 1, 4 1,6 1,8 2,0 0,01 0,1 1 10 HI
[M sun pc -2 ] Galactocentric Radius [R 25 ]
Figure 4: Measured (squares) [START_REF] Bigiel | Extremely inefficient star formation in the outer discs of nearby galaxies[END_REF] and calculated with (39) (solid line) column density of HI for NGC3198.

1) Rotation Curves of spiral galaxies.

2) The peculiar velocities of galaxies in clusters and the problem of the gravitationally bound state of clusters.

3) Gravitational lensing on large structures. As to the first problem, we have shown clearly in present paper, that rotation curves of galaxies can be explained without invoking the hypothesis of the existence of a DM. Actually the RC can be obtained directly from the baryon density distribution if we recognize that the gas movement obeys the diffusion laws.

The second problem deserves more detailed consideration for the following reasons: 2a) In the worst case of large cluster, characterized by size 3 M pc and the velocities dispersion 1000 km/s , we obtain for a galaxy the cluster crossing time t = 3 • 3 • 10 9 = 10 10 yrs that is comparable with the cosmological time. This mean that the cluster does not need to be bounded in the strict sense of the word. [17], [START_REF] De Blok | High-resolution rotation curves and galaxy mass models from[END_REF] and calculated with (22a) (dashed line) rotation curve for NGC7331. Wind tail that corresponds to the HI distribution (see fig. 3) is shown by the horizontal bold solid line. 2b) Strictly speaking the virial theorem can not be directly applied to the cluster if some of its members have masses comparable (and this is the case) with the cluster mass.

2c) The geometry of our universe does not have to be (pseudo-) Riemannian. Moreover, there are serious reasons to believe that we live on the Finsler manifold (see discussion below) and in general case the dynamics for large mass and distances will distinct of the Newtonian one.

All this clearly indicates the weakness of the second argument.

If we talk about the third argument in favor of the DM, first of all, we must remember again that we do not know the metric tensor of the variety on which we live. For this reason it is not entirely correct to use the expression for lensing, obtained for a particular case of Riemannian geometry.

However, let us return to cosmology. It is well known that DM plays a key role in the formation of the observable structure of the universe. Within the framework of the (pseudo-) Riemannian geometry, cosmological time is not sufficient for the observable structure be formed in the absence of the DM. However, now there are several arguments against the need for DM. On the one hand, as it was shown by P. Kroupa et al [START_REF] Kroupa | The failures of the standard model of cosmology require a new paradigm[END_REF], cosmological models based on warm or cold DM are not able to explain observed regularities in the properties of dwarf galaxies. On the other hand, as we have shown, we do not need DM to explain rotation curves of galaxies. All this clearly indicates that the [17,[START_REF] Begeman | HI rotation curves of spiral galaxies[END_REF][START_REF] De Blok | High-resolution rotation curves and galaxy mass models from[END_REF] and calculated with (22a) (dashed line) rotation curve for NGC3198. Wind tail that corresponds to the HI distribution (see fig. 4) is shown by the horizontal bold solid line.

paradigm should be revised. The only reasonable extension of the existing paradigm, which satisfies the principle of the Occam's razor, is the extension of the (pseudo-) Riemannian geometry to the Finslerian one, that gives the necessary time we need to form the observed structure. As it was said by Shiing-Shen Chern, "Finsler geometry is just Riemannian geometry without the quadratic restriction" [START_REF] Chern | Finsler geometry is just Riemannian geometry without the quadratic restriction[END_REF]. Actually the (pseudo-) Riemannian geometry is a very special case of the Finslerian one, and there is no reason for this particular restriction. Moreover, only within the framework of the Finslerian geometry, the cosmological constant appears in a natural way from geometry itself, it has natural explanation and it becomes possible to unify quantum theory and gravity [START_REF] Lipovka | Physics on the Adiabatically Changed Finslerian Manifold and Cosmology[END_REF], [START_REF] Lipovka | Planck Constant as Adiabatic Invariant Characterized by Hubble's and Cosmological Constants[END_REF].

Besides that on the Finslerian manifold the Planck constant calculated from the first principles (with measured cosmological parameters) coincides with it's experimental value up to second digit [START_REF] Lipovka | Physics on the Adiabatically Changed Finslerian Manifold and Cosmology[END_REF] , whereas if it is calculated in the (pseudo-) Riemannian world, we find that the Planck constant differs by factor 3/2 from it's exact value [START_REF] Lipovka | Variation of the fine structure constant[END_REF]. This is more than a serious argument in favor of the Finsler geometry.

The main results of the paper can be summarized as follows: 1) We obtain an analytic expression for the density function of spherically distributed matter, satisfying the following condition: V ⊥ = const and compared it with those found for DM from numerical simulations of rotation curves produced by DM.

2) A new density distribution function is proposed, which allow integrate analytically.

3) On the basis of new factorized density function, which excellently fits the Miyamoto-Nagai one, we obtain a general expression for the galactic RC and its baryon mass.

4) We argue that RCs are formed of two parts. One (inner) is formed by collisionless ideal gas, consisting of stars, and the other (outer) by the real gas, whose motion obeys not only the gravity, but also the laws of kinetics (diffusion equations).

5) The direct and exact relationship between RC and outer gas density function is shown. We calculate the HI column density functions, that correspond to the measured rotation curves for two edge-on spiral galaxies: NGC7331 and NGC3198. The calculated column density are in excellent agreement with the observed.

6) The total mass of three spiral galaxies is calculated. Our evaluation for the Galaxy is M G = 23 • 10 10 M ⊙ . For NGC7331 we find M 7331 = 37.6 • 10 10 M ⊙ and for NGC3198 the total mass is M 3198 = 7.7 • 10 10 M ⊙ .

  Figure 1 suggests this comparison for the plane z = 0. Parameters for the best fitting of the M-N function by the expression (18) are as follows: γ = 30, α 1 = 0.024 , β 1 = 12 , α 2 = 3.7 , and β 2 = 10000.

Figure 1 :

 1 Figure 1: Baryon density as a function of distance R for bulge + disk model of the Galaxy. The density of Miyamoto -Nagai (solid line), and density given by the expression (18) (dashed curve) are compared.

Figure 2 :

 2 Figure 2: Rotation curve obtained with density (18) for the Galaxy.

Figure 5 :

 5 Figure5: Measured (squares)[17],[START_REF] De Blok | High-resolution rotation curves and galaxy mass models from[END_REF] and calculated with (22a) (dashed line) rotation curve for NGC7331. Wind tail that corresponds to the HI distribution (see fig.3) is shown by the horizontal bold solid line.

Figure 6 :

 6 Figure 6: Measured (squares)[17,[START_REF] Begeman | HI rotation curves of spiral galaxies[END_REF][START_REF] De Blok | High-resolution rotation curves and galaxy mass models from[END_REF] and calculated with (22a) (dashed line) rotation curve for NGC3198. Wind tail that corresponds to the HI distribution (see fig.4) is shown by the horizontal bold solid line.

  (is shown by squares with error bars).

	The parameters of the model are α * 1 = 0.2317 , β * 1 = 0.112 , α * 2 = 6.358 , β * 2 = 28.8, α * 3 = 7.005 , and β * 3 = 1440.
		300
	Rotation curve (Km/s)	100 150 200 250
		50
		0,1	1	10	100
			Distance R (kpc)

Acknowledgments

I would like to express my gratitude to Dr. I.V. Gosachinsky for useful discussions.

Bibliography