Application and experimental validation of new computational models of friction forces and rolling resistance

In the process of modelling of the Celtic stone rotating and rolling on a plane surface, different versions of simplified models of contact forces between two contacting bodies are used. The considered contact models take into account coupled dry friction force and torque, as well as rolling resistance. They were developed using Padé approximations, their modifications, and polynomial functions. Before the use of these models, some kinds of regularisations have been made, allowing to avoid singularities in the differential equations. The models were tested both numerically and experimentally, giving some practical guesses of the most essential elements of contact modelling in the Celtic stone numerical simulations. Since the tested contact models do not require the space discretisation, they can find application in relatively fast numerical simulations of rigid bodies with friction contacts.

of admissible loadings of the contact for which the developed sliding does not occur. For determination of the sliding direction, they used the properties of the convex zone of admissible loadings of the contact, i.e. the fact that its normal direction determines the components of sliding [9]. But in this case, the approximation by the use of the ellipsoid leads to much higher errors.
Zhuravlev [30] developed essentially the results of Contensou for circular contact area with Coulomb friction law and fully developed sliding, by giving exact analytical expressions for friction force and moment. Since the obtained analytical results are not convenient in practical use, he also proposed the corresponding firstorder Padé approximants as effective approximation of the integral model of friction [31]. For the same contact, the coupled friction model was approximated by the use of Taylor expansion of the velocity pseudo potential [21]. The same group of authors proposed for this contact an ellipse as an alternative efficient approximation of the relation between the magnitude of friction force and moment [27]. The piece-wise linear approximation of integral friction model for elliptic contact area and the Hertz stress distribution was presented in paper [15]. Kireenkov [15] constructed the coupled model of friction forces and rolling resistance for circular contact area with Coulomb friction law and fully developed sliding. He used the second-order Padé approximants, which are more accurate and suitable for qualitative analysis. The rolling resistance was modelled as a result of distortion of initial Hertzian stress distribution.
In the work [19], the authors presented some generalisations and modifications of Padé approximants and applied them in modelling of friction forces for fully developed sliding on circular contact area, with constant contact pressure. The presented new family of models can be understood as a bridge between Padé approximations and models proposed by Möller et al. [27]. The continuation of the project [19] is presented in work [20], where an integral model of dry friction components is built for a general shape of plane contact area. Then, the special model of stress distribution over the elliptic contact patch is developed, being a generalisation of Hertzian distribution, with addition of special distortion related to the resistance of rolling. At the end, some approximations of the resultant friction forces are exhibited, being piece-wise polynomial functions as well as further generalisations and developments of Padé approximants.
A good object for testing coupled mathematical models of dry friction components and rolling resistance for elliptic contact area is the Celtic stone, also known as wobblestone or rattleback. It is usually a semi-ellipsoidal solid (or another kind of body with a smoothly curved oblong lower surface) with a special mass distribution. Most celts lie on a flat horizontal surface and when set in rotational motion around the vertical axis can rotate in one direction only. The imposition of an initial spin in the opposite direction leads to transverse wobbling and then to spinning in the preferred direction.
The Celtic stone with its special dynamical properties was an object of investigation of many researchers, and the first scientific publication on this subject appeared at the end of the nineteenth century [28]. Dissipationfree rolling without slip is one of the widely used assumptions in modelling of the celt [4,5,23,28]. Walker [28] pointed out the non-coincidence of the principal axes of inertia and the principal directions of curvature at the equilibrium contact point as essential in explanation of the wobblestone properties. An analysis of the linearised equations of the model assuming continuous slipping (quasi-viscous relation between the friction force and velocity of the contact point) is presented by Magnus [24]. Another model which takes into account dissipation, but is far from reality, is analysed using the asymptotic perturbation theory [6]. Some authors proposed a model, which assumes rolling without slip and viscous damping (torque about all three axes) [13]. Garcia and Hubbard [8] presented a more realistic modelling with aerodynamic dissipation and slip with dry friction force. They completed the model with an experimental validation. Markeev [25] performed the perturbation analysis of local dynamics around the equilibrium points of the model assuming the absence of friction. He also gave the experimental verification. The modelling of celt closer to reality is proposed by Zhuravlev [32], who assumed the possibility of slip, but in contrast to all the earlier works, applied the spatial Coulomb-Contensou friction model with linear Padé approximations for circular contact patch [14,30,31]. However, since the friction force is the only way of dissipation in the proposed model, the time of the wobblestone motion (until rest) is unrealistically long. In the papers [16,17] there is presented an attempt to more realistic modelling of the Celtic stone dynamics, where both friction force and torque are taken into account, as well as rolling resistance. The Coulomb-Contensou model and the first-order Padé approximants are applied to model the frictional coupling between translational and rotational motion of the circular contact area. The dependence of the dynamical response on initial conditions and the system parameters is investigated, showing essential differences in the properties of the proposed model, when compared with the previous simulation models of the celt. The complete set of hyperbolic tangent functions approximating friction components coupled with rolling friction for circular contact area is proposed and then applied in the wobblestone dynamics modelling in the work [18]. In the works [1,2], the series of the Celtic stone models, with different levels of complexity of the model of the contact forces for the elliptic contact zone (from the independent friction force and torque components, to the second-order Padé approximants), are investigated. The simulations are compared with the corresponding results obtained by the use of the model, where the friction components of the integral model are approximated by more precise, but much more complicated approximations based on the piece-wise polynomials. The untypically forced wobblestone, i.e. situated on a harmonically vibrating base, is investigated numerically in the work [3], where rich bifurcation dynamics is detected and presented. Section 2 of the present work is devoted to mathematical modelling of the Celtic stone as a rigid body with continuous contact with the plane surface. Section 3 presents the corresponding approximate coupled models of dry friction components (spatial force and torque), rolling resistance and their implementation in the wobblestone model. Sections 3.1 and 3.2 are in fact a kind of an extract from the paper [20]. In Sect. 4 we present results of numerical and experimental investigations: parameter estimation results, analysis of different versions of the contact model, comparisons between simulation and experimental solutions. Section 5 gives some final remarks.

Differential and algebraic equations of motion
The wobblestone as a semi-ellipsoid rigid body with the geometry centre O and mass centre C (with the relative position defined by vector k) touching a rigid, flat and immovable horizontal surface π (parallel to the X 1 X 2 plane of the global immovable coordinate system G X 1 X 2 X 3 ) at point A is presented in Fig. 1 [1][2][3][16][17][18].
The differential equations of motion can be written as follows: Bd where m is the mass of the celt, B is the tensor of inertia of the solid in mass centre C, v is the absolute velocity of mass centre C, ω is the absolute angular velocity of the body, r C = r GC is the vector defining the position of mass centre C with respect to the origin G of the global coordinate system,N =N n is the normal reaction of the horizontal plane, n is the unit vector normal to the plane X 1 X 2 ,T Fig. 1) is the sliding friction force in the point of contact A,M Fig. 1) are the dry friction and rolling  (1) represents the law of motion of the mass centre of the body in the absolute coordinate system, while Eq. (3) exhibits the law of variation of the angular momentum expressed in the local frame, where the notationdu/dt stands for the time derivative of vector u in the movable body-fixed coordinate system. Equation (4) represents the relations between the components of ω in the local reference frame and derivatives of Cardan angles describing the orientation of the body by the following sequence of three rotations about the axes of the local coordinate system C x 1 x 2 x 3 (of axes parallel to the geometric axes of the ellipsoid): x 3 (by an angle ψ), x 1 (by an angle θ ) and x 2 (by an angle ϕ).
For any vector u = u x 1 e x 1 +u x 2 e x 2 +u x 3 e x 3 , where e x 1 e x 2 and e x 3 are the unit vectors of the corresponding axes of the local coordinate system C x 1 x 2 x 3 , we use notations u x i = u i and e x i = e i (i = 1, 2, 3). Moreover, we denote the components of inertia tensor B in the local reference frame as follows: where B i = B x i is the moment of inertia of the body with respect to axis C x i (i = 1, 2, 3), while B i j = B x i x j is the corresponding centrifugal moment of inertia (i, j = 1, 2, 3 and i = j).
Since the numerical simulations are performed in both the local and global reference frames, one needs the following rule of transformation between the corresponding coordinates of a vector u: where Applying the above transformation to the coordinates of vector n, we get In order to find the relation between vectors r and n, one can use the following ellipsoid equation: (where a 1 , a 2 and a 3 are the semi-axes of the ellipsoid along the axes x 1 , x 2 and x 3 , respectively) and the condition of tangent contact between the ellipsoid and horizontal plane where δ < 0. The solution to Eqs. (9) and (10) with respect to the components of r and δ in the C x 1 x 2 x 3 coordinate system reads where The set of Eqs. (1)-(4) consists of 12 scalar first-order ordinary differential equations with 13 unknown functions of time: v X 1 , v X 2 , v X 3 , r C X 1 , r C X 2 , r C X 3 , ω 1 , ω 2 , ω 3 , ψ, θ , ϕ andN . The missing equation is the following algebraic expression: which follows from the fact that the point A always lies in the plane π. Equations (1)- (4) and (13) form now the set of differential-algebraic equations of index 3. Differentiating the condition (13) with respect to time in the global coordinate system G X 1 X 2 X 3 , one gets Since dr C /dt = v,dk/dt = 0 and (dr/dt) · n = 0, Eq. (14) takes the following form: which expresses the fact that velocity v A = v + ω × (r + k) lies in the plane π. Equations (1)- (4) and (16) are now the differential-algebraic equations of index 2.
Differentiating the algebraic condition again, one obtains Since we finally get the following differential equation: Since in the above expressions the local derivatives of vector r appear, one need to differentiate the expressions (11) with respect to time and then do the same with the local components (8) of vector n. Now we have the set (1)-(4) and (19) of 13 scalar differential-algebraic equations of index 1, with 13 unknown functions of time, but with only 12 derivatives of unknown functions. One can solve them with respect to the corresponding derivatives and the normal reactionN and obtain the set of 12 differential equations with 12 unknown functions. During numerical simulations of the Celtic stone, one should be very careful when choosing the numerical method and initial conditions, since, especially during the long time simulations, there is a risk of violation of two algebraic conditions (13) and (16).
Note that we have chosen the local coordinate system with the axes parallel to the geometric axes of the ellipsoid, instead of central and principal axes of inertia. In the second case the tensor of inertia takes a convenient form of diagonal matrix. However this advantage is lost since the unknownN appears in the right-hand side of Eq. (3) (note that functionN appears also in the expressions forT . Moreover, the relation between the vectors r and n is simpler in the chosen coordinate system C x 1 x 2 x 3 .

Principal curvatures of the ellipsoid
For some developed later contact models we will take into account time variation of the curvatures of the wobblestone surface in the point of contact with the base. While limiting considerations to the surface of revolution, the following parametric description of the wobblestone surface can be proposed: where for parameters 0 ≤ p ≤ π and 0 ≤ q < 2π. Now the principal curvatures can be expressed as follows [10]: where κ 1 = R −1 1 is the curvature in the meridional direction (in the section including the Ox 1 axis), whereas κ 2 = R −1 2 is the curvature in the corresponding transverse direction (R 1 and R 2 are the corresponding curvature radii).
From Eqs. (22) and (21) we get Equations (20)-(21) give the relation r 2 which leads to the following form of Eqs. (23): The curvature ratio of ellipsoid κ = κ 2 /κ 1 reads where κ > 1 for a 1 > a 2 and −a 1 < r 1 < a 1 . In the case of the contact between the Celtic stone with the plane surface, the principal curvatures and the corresponding principal directions of the ellipsoid surface are also the relative principal curvature and relative principal directions of the contact. As mentioned in the Sect. 1, the non-coincidence of the principal directions of curvature and the principal axes of inertia at the equilibrium contact point is responsible for the typical dynamical properties of the Celtic stone. This can be realised in the two ways: (i) an asymmetrical base of the stone with the skewed principal directions of curvature and (ii) a symmetrical base, but special mass distribution is obtained by the use of additional elements mounted on the stone. Since we assume a geometry of the body in the form of an ellipsoid, the wobblestone investigated in this work is of the second kind.

Friction and rolling resistance modelling
3.1 Contact pressure distribution, rolling resistance and integral model of friction forces Figure 2 exhibits the non-dimensional elliptic contact area F with the co-ordinate system Ax y, where the vector n faces the wobblestone. The real co-ordinates of any point arex =âx andŷ =â y, whereâ is the real length of the major semi-axis of the elliptic contact.
Based on the work [20], we assume the following non-dimensional distribution of the contact pressure between the Celtic stone and the base plane: where The function σ 0 ρ represents a certain circularly symmetric stress distribution on the unit circle F , where ρ is a distance from the centre of the area F . We assume here that , whereN is the normal component of resultant real force of interaction between bodies. The symbol d a can be interpreted as a non-dimensional rolling resistance coefficient along the major axis direction, while the parameter b r describes the corresponding rolling resistance coefficient b r d a along the minor axis of the elliptic contact area. Since we assume that a relatively rigid wobblestone rolls on the deformable table, the quantities v r x and v r y are the components of the non-dimensional velocity v r =v r /â of relative motion of the point A over the plane π (wherev r is its real counterpart). In order to avoid the singularities and the set-valued force laws, we have added here a small parameter ε R under the square root in Eq. (29). The proposed model (27) of the contact pressure distribution should be treated rather as artificial mathematical construction, being a kind of generalisation of the Hertzian stress distribution on an elliptical contact area, with additional distortion related to the rolling resistance.
The non-dimensional rolling resistance, additionally denoted as "A model", results from the position of the centre of the contact pressure distribution and reads where The rolling resistance moment, which appears in Eqs. (1)-(4), readŝ where X stands for a version of the model. Further versions of the rolling resistance models are defined in Sect. 3.4. When the shape of the contact area assumed in the model is close to the real one, and the material is isotropic, we expect that b r is close to 1. It is consistent with the hypothesis of the linear hysteresis, where the real rolling resistance coefficient is proportional to the length of the contact. Let us also note that we do not take into account the relative angular motion of the deformation zone and its influence on stress distribution and rolling resistance.
Let us assume the fully developed sliding on the contact area shown in Fig. 2. Let the relative motion of the Celtic stone surface in the contact zone F with respect to the fixed table be a locally plane motion. Then the motion can be described by the use of the dimensionless linear sliding velocity v s =v s /â (wherev s is its real counterpart) at the centre A of the contact and angular sliding velocity ω s . Let us also assume the classical Coulomb friction law on each element d F of the contact zone F, with the non-dimensional elementary friction force dT s = dT s /(μN ) acting on the celt, and its moment dM s = dM s /(â μN ) with respect to the pole A, where dT s is the real elementary force, dM s -the corresponding real moment and μ-dry friction coefficient. Assuming the resultant friction force and torque as T s = −T s x e x − T s y e y and M s = −M s n, one gets In Eq. (32) v s denotes the magnitude of the sliding velocity v s , while ω s is the projection of the vector ω s on the vector n.

Approximate models of friction forces
In this section we present a series of different models approximating the integral model of friction (32), for the contact pressure distribution (27), which do not require the space discretisation during the simulation of the wobblestone. All the exhibited models derive from those presented in the work [20], but in order to avoid the singularities and the set-valued friction laws, they are additionally smoothed by means of introduction of a small numerical parameter ε T . In general, we can writê where X i stands for a specific version of the contact model.
The simplest approximation to be tested, denoted as A 0 , is one ignoring the coupling between friction force and torque, where and where v s x and v s y are components of the vector v s = v s x e x + v s y e y along the corresponding axes.

Generalisations of Padé approximants
In the work [20] the following approximation of the integral model (32) was proposed: for f = T sx , T sy , M s , which can be understood as a kind of generalisation of a certain type of the Padé approximants. The coefficients a f,i = a f,i (ϕ s , sgn(ω s )) can be found from the following conditions: for i = 0, ..., n 1 , where f = T sx , T sy , M s .
Here we propose the following regularised version of the approximation (36): for f = T sx , T sy , M s , where the variable v s has been replaced by the quantity v sε , defined by Eq. (35).
The coefficients a f,i = a f,i (ϕ s , sgn(ω s )) of the approximation (36) while for n = 3 and n 1 = 1 we have where one has assumed that m T sx = m T sy = m T s and b T sx = b T sy = b T s . In Eqs. (39)-(40), the parameters and functions defined in the following way have been used: where K (e) and E(e) are the complete elliptic integrals of the first and second kind, while e = √ 1 − b 2 is the eccentricity of the contact. For numerical evaluation of the elliptic integrals during the simulation process see Refs. [20,26].
Assuming the general case of the model of the contact pressure distribution presented in Sect. 3.1, all the parameters (b, I 0 , I 2 , I 3 , m T s , b T s , m M s , b M s , b r and d a ) occurring in the non-dimensional models of friction forces (39)-(40) (see also (29)) will be identified in Sect. 4 from the experimental data.

Piecewise polynomial approximations
In the work [20] the following piecewise polynomial model is constructed: In the present work we will use the following regularised version of the approximation (42): applied to the Hertzian nominal stresses σ 0 ρ = 3/2π 1 − ρ 2 . The coefficients a f ε ,i and b f ε ,i are constructed from the functions a f,i and b f,i given in the work [20], where the expressions cos ϕ s and sin ϕ s are replaced correspondingly by the following quantities: We assume that the parameters u 0, f in the model f are the same as in the approximations f (W ) for Hertzian nominal stresses [20]. Finally, the model f (A H ) ε possesses the same properties as the model f (W ) , in the limit ε T → 0.
In the model A H we will assume, besides the Hertzian nominal stresses, the shape and size of the contact according to the Hertz theory and the approximated model presented in Sect. 3.3, for the actual (changing in time) relative principal curvatures according to the model described in Sect. 2.2. The model A H will be tested in Sect. 4, but we have developed it also in order to substitute the full integral model (since it can be assumed as an almost ideal approximation of the integral friction model for the Hertz nominal stress distribution) in further numerical experiments and investigations of importance of particular elements of the approximate friction models in other cases of the Celtic stone dynamics or other dynamical system with elliptical friction contacts.

Shape and size of the Hertz contact
According to the Hertz theory [12] the shape (eccentricity e or b parameter) of the elliptic contact patch is determined by the following equation: where κ = κ 2 /κ 1 is the ratio of relative principal curvatures κ 2 and κ 1 , where κ 2 > κ 1 and the longer axis of the contact area has the direction of curvature κ 1 . Because Eq. (45) can be solved only numerically, the following approximate solution is given in the work [12]: It turns out that better results can be obtained by the use of approximation b (a 2 ) = κ −0.652 , where we have taken a bit different exponent. Even better results are given by the following approximate expression: where for the assumed form of the function the numerical values of the coefficients have been optimised. Absolute errors of approximations b (a i ) as functions of κ argument along with the exact plot of b (κ) are presented in Fig. 3. The absolute size of the elliptic contact area is determined by the following equation [12]: whereN is the normal load of the contact, while E * is the material constant expressed as follows: where E 1 , E 2 , ν 1 and ν 2 are Young's moduli and Poisson's ratios of the materials of two contacting bodies, respectively. Taking into account thatb 1/2 = (bâ) 1/2 = (1 − e 2 ) 1/4â1/2 , the relationship (48) can be presented in the following form:â where F(e) is the following function: Although the above expression is undetermined for e = 0, the limit lim e→0 F (e) = 1 exists. One can assume some threshold value of eccentricity, e.g. e = 10 −3 , for which F 10 −3 ≈ 1 + 2, 5 · 10 −7 and below which the value of function F(e) can be assumed as 1. Function F(e) is ascending on the interval 0 < e < 1, and for example F 2 · 10 −4 ≈ 1 + 10 −8 . Computing numerically the values of function F(e) from formula (51) for smaller arguments, the numerical instabilities can be encountered. The results (47) and (50) are used during the numerical simulations presented in Sect. 4, as elements of the contact model A H , approximating shape (b) and size (â) of the elliptic contact area.

Implementation of the contact forces models
In order to apply the developed models of the resultant friction forces and rolling resistance in the numerical simulations of the wobblestone mathematically modelled in Sect. 2, we need some relations between the corresponding coordinate systems and vectors.
The versors of the coordinate system fixed to the contact patch are defined in the following way: where is the projection of versor e 1 of the axis x 1 onto the plane π (cf. Fig. 1). Here, we make an assumption that the axis x 1 is never perpendicular to plane π. The non-dimensional sliding, spinning and rolling velocities are given by the following relations: larger. Three small steel elements 5, isolated from the celt thermally and playing the role of measurement points, were placed on the stone. They were heated up before the experiment and then observed by a thermovision video camera. Point P 3 was located on top of rod 4 attached to the stone in such a way that it coincided with the x 3 axis of the local reference frame exhibited in Fig. 1.
The measurements were performed using a single thermovision video camera of 100 Hz frequency and 120 × 160 definition, located at a certain height above the stone with the axis of lens perpendicular to the plane base. The absolute position of the stone was determined uniquely (assuming the permanent contact between the stone and the base) by the coordinates of three points P 1 , P 2 and P 3 observed in the camera image. The coordinates of the measurement points in the local body-fixed reference frame were r O P 1 1 = −97 ± 0.5 mm, r O P 1 2 = 0 ± 0.5 mm, r O P 1 3 = 4 ± 0.5 mm, r O P 2 1 = 97 ± 0.5 mm, r O P 2 2 = 0 ± 0.5 mm, r O P 2 3 = 4 ± 0.5 mm, r O P 3 1 = 0 ± 0.5 mm, r O P 3 2 = 0 ± 0.5 mm, r O P 3 3 = 151 ± 0.5 mm. Figure 5 shows an exemplary frame from the experimental movie compared to the corresponding one obtained from animation based on the numerical simulation.
In order to compare the positions of points P 1 , P 2 and P 3 obtained from numerical simulations to the corresponding experimental locations, we express the corresponding absolute positions of the points in the following way: Then they are projected centrally, with the centre at the focal point of the camera, onto the plane π.
During the process of estimation of the parameters and initial conditions, the following objective function: describing quantitatively the fitting of numerical simulation to the real system behaviour will be minimised.
In the above formula X for the corresponding coordinates along the axis X k of the projections P j andP j of the measurement points obtained in the simulation process (P j ) and experimentally (P j ), in the r -th series (a single system solution with one set of initial conditions) at time where t is the sampling period. Moreover, the symbol M denotes the number of different solutions compared. The process of minimisation of the objective function (58) is a highly dimensional (parameters of the system and initial condition for each series are estimated) nonlinear optimisation problem which we solve using the simplex method with certain modifications proposed in the work [29] and minimising the problem of local minima by bootstrap restarting. Let us assume an objective function in the form F o (p, y), where y is the vector of experimental data and p is the vector of parameters (including initial conditions). The bootstrap restarting method can be then presented as follows [29]: 0. Using the starting vector p 0 , find a local minimum of F o (p, y):p 0 . Then repeat the steps 1-3 for i = 1, . . . , K b . 1. Create bootstrap resample y * i of original data y. Using the starting vectorp i−1 , find a local minimum of F o p, y * i : p * i . 2. Using the starting vector p * i , find a local minimum of F o (p, y): The simplex method is used to find a local minimum. We obtain the bootstrap resample of the same size as the original dataset using random sampling with replacement from the original dataset. We choose the number of necessary iterations K b by observing the values of the objective function F o p i , y against i. If no significant changes of the objective function take place, we can assume that the global minimum is reached with relatively high probability and accuracy.
In the experimental investigations we used a rubber plate as a base for the celt. Therefore, we can assume that the stone was absolutely rigid in comparison with the base plane (the contour resistance for the celt can be neglected). We used three different experimental solutions in the identification processes for different versions of the contact model. The final fitting of three selected models A 0 , A 2 and E 0 to the experimental data is presented in Fig. 6 exhibiting Euler angles for three different solutions, respectively. Table 1 exhibits parameters of the selected models (the parameters used in the regularisation process are ε T = 10 −2.5 and ε R = 10 −1 ). The parameters with the values typed in italics are not independent in a given model (they are assumed arbitrarily and do not participate in the identification process). Figure 7a shows corresponding final values of the objective function obtained in the minimisation process for the experimental data. One can note that models A 1 , A 2 , B 1 , B 2 , C 1 , C 2 , D 1 and D 2 cover the dynamics of the celt dynamics almost in the same degree (the smallest value of the objective was obtained for model A 2 : F O = 1.70 mm 2 ), while other ones lead to significantly worse results. It means that such properties like orthotropy of the rolling resistance as well as coupling between friction force and torque play an essential role in the modelling process of the Celtic stone, while other elements are less important. Figure 8 presents other quantities obtained from the simulation of model A 2 (component ω 3 of the angular velocity, normal reaction of baseN and total mechanical energy E m ), for the two first solutions only. Finally,  Fig. 9 illustrates the friction force components and friction torque for models A 0 and A 2 , where one can observe more visible differences in friction torque.
The identification process was also performed for model A H representing Hertz stress distribution over the contact area. The obtained value of objective function F o = 2.30 mm 2 is larger than the value corresponding to model A 2 since assuming the Hertz contact we have a smaller number of the model parameters. We developed this model in order to make further numerical experiments and investigations of importance of particular elements of the friction model in other cases of the Celtic stone dynamics and other dynamical systems with elliptical friction contact (since the developed piecewise polynomial approximation can be assumed as an  almost ideal approximation of the integral friction model for the Hertz normal stress distribution). We also performed a numerical experiment of fitting the celt models with different approximant friction models to the three numerical simulations obtained by the use of model A H (the same simulations as these obtained in the identification from experimental data). The initial conditions in this case are known from numerical simulation, and they are not identified (it is a reason of larger errors of the model E 0 ). The results are presented in Fig. 7b, and we can observe that they are qualitatively the same as the result of identification from the experimental Fig. 9 Components of friction models A 0 and A 2 for the simulation results corresponding to the first two experimental solutions (for the legend see Fig. 6) data (Fig. 7a). Models A 1 , A 2 , B 1 , B 2 , C 1 , C 2 , D 1 and D 2 lead to almost the same results, and the smallest value of the objective function obtained for model A 2 is F o = 0.0048 mm 2 . We can conclude that the model D 1 (as the simplest from the group A 1 , A 2 , B 1 , B 2 , C 1 , C 2 , D 1 and D 2 ) is sufficient for modelling of the friction model in the celt if the real stress distribution is close to Hertz distribution and friction laws are close to Coulomb friction law for each element of the contact area.

Concluding remarks
The present work makes use of the approximate models of coupled friction and rolling resistance in the case of elliptic contact, developed in the work [20]. The corresponding models are applied in the Celtic stone modelling and then tested both numerically and experimentally. Both the proposed model of the celt and its simulations are very realistic, when compared to the majority of earlier works on it: they take into account correct models for linear and torsional friction as well as rolling resistance. In particular, the present work shows the importance of the following elements of the modelling of typical dynamics of the celt: coupling between friction force and torque and orthotropy of rolling resistance related to the elliptic shape of contact. The assumption of the elliptic shape of the contact turns out to be less important.