Rolling resistance modelling in the Celtic stone dynamics

In the work there is presented a mathematical model of a rattleback lying on a horizontal plane, with special attention paid to modelling of the contact forces. Friction force modelling concerns a certain class of problems, where it is possible to assume fully developed sliding on the contact zone. The models are based on the integral expressions assuming the classical Coulomb friction law valid at each point of the contact. In order to obtain a higher simulation speed, a special class of approximations if integral functions is used, being some kind of generalization of Padé expansion. The rolling resistance is modelled as a result of special distortion of contact pressure distribution and takes into account the rolling direction and the elliptic shape of the contact. The work uses previous results of the authors, but more attention is paid to the rolling resistance model, which is extended and analysed in more detail. The models are validated experimentally and their various elements are tested.


Introduction
Modelling and numerical simulation of mechanical systems with frictional contacts in their general case require discretization of the deformable contacting bodies near the contact zone or at least discretization of the contact surface [1][2][3]. For the rolling contact one can use the well-known CONTACT software [2]. Space discretization and use of finite element methods lead, however, to high computational costs and long times of numerical simulations that are not always acceptable. It is the reason of the development of simplified models and algorithms [3][4][5].
There are special situations, where the local deformations of the bodies are negligible in comparison with their relative velocities in the contact area, and one can assume a fully developed sliding over the contact. In the case of planar contact one can assume relative motion of the contacting bodies in the contact zone, equivalent to the planar motion of a rigid body. Modelling of the resulting friction force and moment then leads to integral expressions over the contact area [6]. In order to obtain higher speeds of numerical simulations, scientists still searched for reduced and simplified models. Zhuravlev [7] applied Padé approximants in the construction of expressions for friction force and moment in the case of fully developed sliding over a circular contact area with Hertzian contact stress distribution and Coulomb friction law. These results were then developed and the models of mutually coupled friction forces and rolling resistance were obtained [8], where initially Hertzian contact stress distribution is distorted in a special way in order to take into account the rolling resistance [9]. An extension of this approach is presented in the work [10], where the generalized forms of approximations of the resulting friction force and moment as well as rolling resistance for elliptic contact shape are presented. Rolling friction is there modelled as a resistance against motion of the area of deformation over the body [11] and, ignoring the problem of rotation of the elliptical region, is consistent with the results obtained under the assumption of the linear hysteresis [12][13][14][15]. Karapetyan [16] proposed a two-parametric model of friction forces for spherical contact between a ball and planar base being a result of deformations of the contacting bodies. This model joins three different kinds of friction (Coulomb, Contensou and rolling friction) and can be seen as a development of Zhuravlev model.
An example of a dynamical system, in which an essential role may be played by the particular elements of the contact models presented in the work [10], is the rattleback (known also as wobble stone or Celtic stone)-usually a semi-ellipsoidal top with a special mass distribution, moving on the plane base and possessing specific dynamical properties. The first scientific work on rattleback was published at the end of the 19th century [17] and it has been an object of interest till today [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34], both from the point of view of dynamics and of modelling. During modelling and analysis of the Celtic stone, rolling without sliding is often assumed, leading to non-holonomic constraints [17, 21, 22, 24-26, 29, 31]. Some authors assume some kind of dissipation, but far from reality [18][19][20]23]. More realistic modelling of friction is applied in the work [27]. Detailed modelling of contact between the Celtic stone and horizontal base, including rolling resistance and using results of the work [10], is presented in Ref. [34], where also experimental validation of particular parts of the model is performed. The same approach is used in modelling and numerical simulations of a billiard ball rolling and sliding on a table [35,36].
In the present work the same experimental data is used as in Ref. [34], but deeper insight into the rolling resistance model is achieved-an extended expression for the rolling resistance additional elements and aspects are investigated.

Governing equations of the Celtic stone
The presented in this section modelling is based on the work [5]. A wobble stone having the form of semi-ellipsoidal rigid body and lying on a horizontal plane π is presented in Fig. 1, where we use the following notation: GX 1 X 2 X 3 is for the global (immovable) coordinate system with X 1 X 2 axes laying on the horizontal plane π ; n is for the unit vector normal to the surface X 1 X 2 ; O is for the geometry centre of the ellipsoid forming surface of the rattleback; C is for the mass centre of the body; Cx 1 x 2 x 3 is for the local (body-fixed) coordinate system Fig. 1 The wobble stone on a horizontal plane of axes parallel to the geometric axes of the ellipsoid and origin at the mass centre; A is for the contact point; k is for the vector defining the relative position of mass centre; r is for the vector defining the relative position of the contact point;N =Nn is for the normal reaction of the horizontal plane of magnitudeN ;T s is for resultant friction force at the point of contact A; m is for the mass of the Celtic stone;M s is for the friction moment applied to the body;M r is for the rolling resistance moment; mg is for the force of gravity.
The presented system is governed by the following set of differential equations: Bd ω where the following notation is used: du/dt -absolute time derivative (in GX 1 X 2 X 3 coordinate system) of vector u;du/dt is for the time derivative of vector u in the movable (body-fixed) coordinate system; v is for the absolute velocity of the mass centre C; ω is for the absolute angular velocity of the body; T is for the tensor of inertia of the solid at the mass centre C (where B i and B ij denote the corresponding moments of inertia with respect to the axes x i and x j ); r C is for the vector indicating the position of the mass centre C with respect to the immovable origin G of global coordinate system; ψ , θ , ϕ is for the angles describing the orientation of the rattleback defined by the following sequence of rotations: about axis x 3 by angle ψ , about axis x 1 by angle θ and about axis x 2 by angle ϕ, ω 1 , ω 2 , ω 3 is for the components of vector ω along the corresponding axes x 1 , x 2 , and x 3 of the local coordinate system. Equations (1)-(2) define the motion of the mass centre C of the Celtic stone in the global coordinate system. Equation (3) represents the law of variation of the angular momentum expressed in the local reference frame. Equations (4)-(6) represent kinematic differential equations expressing relations between components of angular velocity and derivatives of angles defining orientation of the body. Equation (7) is the twice differentiated algebraic condition (r C + r + k) · n = 0 expressing assumption of permanent contact between the rattleback and plane π . The set of Eqs. (1)-(6) is in fact a set of 13 scalar differentialalgebraic equations of index 1. One can notice here 13 unknown scalar functions of time, but only 12 derivatives of these functions occur-there is no derivative of the functionN(t). Equations (1)-(6) can be solved algebraically with respect to these 13 unknown functions of time. One can obtain then the set of 12 ordinary differential equations with 12 unknown functions of time and additional algebraic expression for normal reactionN(t).
The solution to Eqs. (1)-(6) also requires the relation between the vector n and independent variables where n 1 , n 2 and n 3 denote components of vector n along the corresponding axes x 1 , x 2 and x 3 , as well as the following relation between the vectors n and r: where In Eq. (10) a 1 , a 2 and a 3 denote the semi-axes of the ellipsoid along the corresponding axes x 1 , x 2 and x 3 . Since in Eq. (7) the local derivative of the vector r appears, one also needs the corresponding derivative of Eq. (9),ṙ as well as the differentiated form of Eq. (8). In Eq. (11)u stands for the local derivative of the variable u.
Based on the previous work of the authors [10,34], the following model of the resultant contact forces is proposed:T where μ is for the is coefficient of friction, f is for the rolling resistance coefficient, while b i (i = 1, . . . , 5) and ε i (i = 1, 2) are the remaining constant parameters, v is for the is velocity of a point belonging to the rattleback and being in the contact with plane π , ω s is for the component of the angular velocity ω along the vector n, e x and e y are the unit vectors of the axes of the planar rectangular coordinate system associated to the elliptical The elliptical contact area with the corresponding coordinate system and its construction is presented in Fig. 2, where e 1 denotes a unit vector of the axis Cx 1 . Models (12)-(13) represent special case of approximations developed in [10] and tested in modelling of the Celtic stone in [34]. They are a minimal but sufficient form of more general models, allowing one to simulate the rattleback dynamics quickly and realistically. They are based on the assumption of fully developed sliding on the contact area and small enough local deformations resulting in relative motion in the vicinity of the contact being planar motion of rigid bodies. The resultant friction action is reduced to the friction force and moment acting at the centre of the finite contact area. Their models are based on the integral expressions under the assumption of the classical Coulomb friction law valid at each point of the contact zone. In order to obtain higher simulation speed, a special class of approximations is used, being some kind of generalization of the Padé expansion. In general, friction forces are coupled with rolling resistance via contact pressure distribution. A simplified model of pressure distribution is assumed, i.e. the Hertzian model, which is then modified (distorted) in a special way in order to take into account the rolling resistance. The parameters b i (i = 1, . . . , 3) in models (12)-(13) have no direct physical meaning but depend on shape and size of the contact, as well as on the contact pressure distribution and model of friction. They can be identified based on the experimental data. In the case of model of friction force and moment (12)-(13), a circularly symmetric contact pressure distribution is assumed and there is no coupling with the rolling resistance (14). On the other hand, the rolling resistance model (13) assumes an elliptical shape of the contact patch. This is a simplified special case of the contact models chosen based on the work [34].
The model of the rolling resistance (14) is an extended version of the full form of model presented in [10,34], assuming that only deformations of the table π participate in energy dissipation during rolling. This is approximately true in the experiment presented in Sect. 3, where deformations of the Celtic stone can be neglected. The rolling resistance of the model presented in the previous work depends on the direction of the velocity v A and has main directions along the corresponding axes of the elliptic contact corresponding to extremal values of resistance (more precisely the centre of contact pressure distribution lies on an ellipse) and during rolling along intermediate directions the rule is fulfilled that the rolling resistance moment is perpendicular to the velocity v A . Rotation of the elliptical region is ignored in the model. In the present work the model (14) is extended by addition of the parameter b 4 (in Refs. [10,34] equal to 1), which leads to the situation that the moment of rolling friction is not perpendicular to the velocity v A (see Fig. 3). The rolling resistance modelM r is the result of the shift of the contact pressure distribution centre S from the contact centre A. The point S lies on the ellipse of semi-axes equal to f and f b 4 /b 5 of directions corresponding to the directions of the axes of the elliptical contact area. However, the slendernesses of the two ellipses are not necessary the same. The parameters b 4 and b 5 Fig. 3 Ellipses of travel of the contact pressure centre S and end of the rolling resistance vectorM r along with the vector v A of the contact centre velocity have no direct physical meaning but can be interpreted geometrically in the way presented in Fig. 3.
Application of the contact force models (12)- (14) requires the use of the following relations between the corresponding unit vectors (see Fig. 2): is the unit vector e 1 of the axis x 1 projected onto the plane GX 1 X 2 . The remaining quantities occurring in Eqs. (12)- (14) are calculated in the following way: Rolling resistance modelling in the Celtic stone dynamics The parameters ε i (i = 1, 2) was initially introduced in order to eliminate singularities in the expressions (12)- (14) for vanishing sliding or rolling velocities. Their values was small enough to smooth set-valued friction and rolling models but not to change too significantly their properties. However, in this paper the model of rolling resistance is tested for higher values of the parameter ε 2 .

Numerical simulations and experimental validation
In this work there are used results of experiment presented in work [34]. In Fig. 4 there are presented experimental setup (a) and the Celtic stone used during experimental investigations (b). The rattleback is built by skew attachment of aluminium rod 2 (equipped with additional masses 3a and 3b) to the semi-ellipsoidal body 1 made of aluminium alloy EN AW-2017A. On the stone there are situated three small steel elements 5a-c, isolated thermally from the rest of the body, playing the role of measurement points. The points are heated before the experiment and then observed by the use of a thermovision video camera. The measurement points are situated in such a way (the point 5c is placed at the end of special rod 4) that the position of the Celtic stone can be computed based on their positions on the plane image observed in the camera located above the stone perpendicularly to the base-assuming a permanent contact between the body and the table.
The remaining system parameters and initial conditions are estimated using a combination of the Nelder-Mead method and random perturbations of data set by bootstrap restarting. The minimized objective function F o describes the fitting of the numerical simulation to the experimental results and is defined as average squared deviation between the coordinates of the measurement points' projected centrally onto a plane π (i.e. as viewed in the camera image) obtained from the experiment and numerical simulation. For more details see Ref. [34]. Three different experimental solutions are used in the identification process (see Figs. 5, 6, 7).
Eight versions of the model of the contact forces (final values of the parameters and objective function F o are presented in Table 1    One can observe that the use of full model (12)- (14) is reasonable. The greater value of the parameter ε 2 leads not only to better results but also makes the equations less stiff. An especially important element of the rolling resistance model is its dependence on the direction of rolling: the value of the resistance along the corresponding semi-axes of the elliptical contact as well as during rolling along intermediate directions.

Conclusions
The reported results have been obtained by the use the same experimental data as used in the previous work of the authors [34]. The extended model of the rolling resistance is used and new aspects of the rolling friction in the case of an elliptical contact are investigated. The model of the resultant friction force and moment is used in minimal but sufficient form (i.e. more details do not lead to noticeable improvement of the results) based on the results of Ref. [34].
First of all, it has been shown that the rolling resistance is smaller along the shorter semiaxis of the contact (compare model E with models A-D). It has only been confirmed that the vector of the rolling resistance moment is perpendicular to the direction of motion of the deformation region only in a certain approximation. The choice of an accidental model of direction of rolling resistance, like model D, gives much worse results than the model where the rolling resistance is always perpendicular to v A (model C). But extension of the model by addition of the parameter b 4 to the set of parameters undergoing a process of identification (model B) leads to even better results. Another improvement of the rolling resistance model introduced in the present work, which allows one to obtain better agreement with experiment, is that of treatment of the parameter ε 2 as an element not only playing a role of regularization of the model, but also a quantity undergoing the process of identification and allowed one to take any value (model A). Finally, we have confirmed that an important role in the modelling of the resultant friction force and moment is played by the coupling elements between them, i.e. the dependence of the friction force on the angular sliding velocity as well as the dependence of the friction moment on the linear sliding velocity (compare model F with models A-D).
The final conclusion is that the model here presented of the contact forces is very effective. In spite of its simplicity it allows for realistic and reliable numerical simulations of complex mechanical systems with frictional contacts, i.e. the rattleback having been the object of investigation presented in this work.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.