3D projection of the LuGre friction model adapted to varying normal forces

In this paper, we develop an adaptation of the LuGre friction model so as to allow the development of the friction force and its application in any directions on systems subjected to varying normal forces. This is achieved by projecting a modified LuGre model adapted to varying normal forces in 3D along an arbitrary orthogonal system. Consequently, the direction of the friction force is automatically oriented in the correct direction, thus stick, stick-slip, and slip behavior can be represented in all directions. The projected friction model has the following friction features: stick-slip, presliding displacement, frictional lag, varying break-away force, viscous friction, Stribeck effect, and is adapted to varying normal forces. The equivalence of this projected LuGre model with the modified one is proven analytically. The friction model is then applied to simulate the friction on two mechanical systems. The first system consists of a cube sliding on a plane with a transition from stick to slip due to varying normal forces and with a pulling force oriented in multiples directions of the contact plane. The second one is a more complex system consisting of three turbine blades that uses friction to damp their resonance. The results obtained for both systems are consistent with literature.

Amontons and Coulomb [1]. In modern technologies, friction knowledge is indispensable for safe and energy-saving design; consequently, the study of friction has grown significantly over the past 50 years [2]. Friction is an extremely complex mechanism, which involves microinteractions between the surfaces in contact. Some research was conducted to link the microscopic behavior of the friction with the macroscopic one [3,4]. However, in a multibody system the general behavior of friction at the macroscopic level is generally sufficient. Consequently, only friction models representing the macroscopic effect of the friction are employed. These friction models are used in a wide range of mechanisms but are critical in applications such as high-precision positioning [5], tire-road contact [6], wheel and rail contact [7], soft soil contact [8], etc.
By nature all friction models are empirical and only valid in their specific scope [9]. One of the most simple friction models is the Coulomb's one. This model highlights some key features of friction: friction is always in the opposite direction of the velocity, is proportional to the normal force, and exhibits two distinct regimes, stick and slip. These considerations represent the basics of friction and are used in many advanced friction models. With time, more and more friction features were added to models. Friction features and a different friction model commonly used are presented in [10]. In multibody dynamics a friction model aims to reproduce the macroscopic effect of friction, to involve measurable parameters, and to be computationally effective. Comparison of numerical efficiency of the following friction models is presented in [1]: (smooth) Coulomb, velocity-based, Karnopp, Dahl, LuGre, elastoplastic, stick-slip, and Gonthier models. A review comparing the static and dynamic friction models is presented in [11]. A detailed comparison of computational and physical characteristics of 21 different friction models is presented in [12]. Among these friction models, the LuGre model is one the most commonly used in multibody system and motion control. This model only requires six parameters, is computationally efficient, and is able to reproduce most of friction phenomena. This model is a dynamic model often considered as an extension of the Dahl model (Sect. 2.1). As it is a dynamic model, one differential equation must be introduced to represent stick and slip regimes. The classical LuGre model is commonly used in simulation of friction in transmission, and as many other friction models, the classical LuGre model is not suited to be applied in systems with varying normal force. The stability of friction models subjected to varying normal force is studied in [13]. For the LuGre model, a modified version that is able to withstand a varying normal force is proposed in [14] and presented in Sect. 2.2.
In many papers the friction model is applied on a simple system. It generally consists of a block sliding on a horizontal or inclined plane, which makes it easy to predict the direction of the friction force. In real applications, this direction can be difficult to obtain or unpredictable. To overcome this difficulty, there are some generalized friction models. The generalized Coulomb model is presented in [15], a modification of the Coulomb's law to eliminate numerical discontinuity is proposed in [16], and the regularization of the general kinetic friction (GKF) model is presented in Sect. 4. However, as these models are velocity dependent, pure stick cannot be represented. For friction model that can represent pure stick such as the LuGre model, a projection of the equations is needed. The 2D projection of the LuGre friction model in the case of tire-road contact is presented in [17,18]. In this application the 2D projection is in the contact plane (the road), and projection is carried out to take into account the difference in behavior of the tire along the longitudinal or lateral axe. Parameter identification and stabilization methods are presented for nonsmooth multibody systems based on the 2D LuGre friction model in [19].
In this paper, we propose an adaptation of the LuGre model that can be used in multiple directions and on systems subjected to varying normal force. To manage the varying normal forces, we use the modified LuGre model presented in [14] with management of the direction of the slip velocity inspired by [17]. Thus in this paper, we present a set of "projected" equations of the modified LuGre model. This approach allows us to compute the friction force direction by computing each component of the friction force vector while being suited to the case of varying normal forces.

Classical LuGre
At the microscopic level, the contact between two irregular surfaces takes place at several asperities. Dahl [20] experimentally showed that these junctions behave like a spring and proposed to represent the contact as two rigid bodies interacting through elastic bristles. Under a tangential load, these bristles deflect, and the mean deformation of these bristles is the state variable Z (Fig. 1).
The classical LuGre friction model is an extension of the Dahl friction model with an arbitrary steady-state friction characteristic that includes the Stribeck effect. Thus it describes the stick-slip effects. In comparison with Dahl's model, the LuGre one represents the junction between two bodies as a spring-damper system and includes a viscous effect [21]. It can represent the following features: stick-slip, presliding displacement, frictional lag, varying break-away force, viscous friction, and the Stribeck effect. The classical LuGre model is a dynamic model where the friction force (F ) is computed as The function G(V ) depends on many factors such as material properties, temperature, lubrication, etc. This function does not have to be symmetrical, and thus a direction-dependent behavior can be modeled. One classical form of this function is  [23], Dahl [20], Bliman and Sorine [24], Leuven [25], and generalized Maxwell slip [26] LuGre GKF Karnopp and When N is the normal force, F k and F s are the dynamic and static Coulomb friction forces, and μ k and μ s are their respective friction coefficients, V st is in relation with the Stribeck velocity, and α ranges from 0.5 to 2 [22].
The steady-state behavior is governed by the following expressions: We can observe that Eq. (8) is equivalent to the general kinetic friction (GKF) model (Sect. 4).
A comparison of the features of the classical LuGre friction model with other commonly used friction models, namely GKF [5], Karnopp [23], Dahl [20], Bliman and Sorine [24], Leuven [25] and generalized Maxwell slip [26], is presented in Table 1. The LuGre model is often used in control applications because it is a relatively simple model that captures the essential properties of friction. It only requires six parameters that can be fitted by experiments.
The limitation of the classical LuGre model occurs when the normal force varies in function of time. Two problems are observed: 1. When the system is under a tangential load and the velocity is null, changing the normal force results in changing the values of F k and F s , so changing the value of G(V ). However, as V = 0, the state variable Z is not influenced, resulting in an unchanged friction force and an incorrect behavior.
2. When the normal force tends to zero, the friction forces (F k and F s ) tend to zero too. As these values are used in the computation of G(V ), it also tends to zero. As G(V ) appears in the denominator of the computation ofŻ, we can observe some numerical trouble.

1D modified LuGre: varying normal force
To represent friction in a system subjected to varying normal force, a modified LuGre model was developed. It was first introduced in applications such as tire-road contact [6] followed by some mechanisms [27].
In [14] the modified version of the LuGre model is defined aṡ where σ M 0 is a constant parameter, which can be interpreted as the aggregate stiffness per unit of normal force [14], and σ M 1 and σ M 2 are constants. In the steady state the equations are The main difference with the classical model comes from Eq. (10). Instead of using the Coulomb dynamic (F k ) or static (F s ) forces, their corresponding friction coefficients μ s and μ k are used. In Eq. (11), the friction force is directly proportional by an "instantaneous" friction coefficient to the normal force. The state variable Z is now independent of N . This way of writing the LuGre equations resolves the problems described in Sect. 2.1, but in the case of varying normal forces during the stick regime, it can introduce some drift [14]. In the case of a constant normal force N equal to N E , the classical and modified models are equivalent if and only if [14]

Projected modified LuGre
In multibody dynamics software, the computation of the friction force can be difficult to perform. If the direction of the friction force is unknown in pure stick as there is no relative velocity between the surfaces, then the friction force cannot be correctly oriented. Also, the specificity of a system can lead to unpredictable direction of the friction force. Furthermore, multiple contact points are sometimes needed to allow rotation of the bodies, and in this case the orientation of the friction force can be very difficult to predict. The LuGre equations presented previously (Sects. 2.1 and 2.2) are not usable in these conditions because they imply knowledge of the direction of the friction force. In this paper, we propose to decompose the modified LuGre equations along any orthogonal coordinate system. The modeling of a contact between surfaces through a contact between points attached to one of the surfaces and the other surface is widespread (for example, in finite element  . For illustration, let us consider a contact that takes place between a point A attached to body j and a plane attached to body i (Fig. 2). The normal of the plane is denoted by n, r A is the coordinate vector of point A with respect to frame j , r P is the coordinate vector of point P with respect to frame i, and the interbody penetration δ is computed as so that the penetration rate can be retrieved from the relative velocity V computed at middle of the penetration zone (point M on Fig. 3): The relative velocity is decomposed along the tangential ( V t ) and normal ( V n ) components, The normal force N is computed at each contact point according to the Hunt-Crossley formula [28] N = K contact δ P K + D contact δ P D dδ dt (19) with • the contact stiffness computed as where A is the area of contact, and σ is set to 3 · 10 −9 m 2 N 1/2 [29], • contact damping D contact set to 4.88 · 10 5 N s/m 3/2 [29], • interbody penetration δ, • fitting exponential coefficients P k = 2 and P D = (P k − 1)/2 [29].
This can be achieved in any orthogonal system, it ensures to have the sliding velocity oriented parallelly to the contact plane. The projection of the modified LuGre equations along an arbitrary orthogonal coordinate system allows us to compute the three components of the friction force (F x , F y , and F z ), so that the resulting force F is correctly oriented. The indices indicate the projection along the corresponding axes (x, y, or z). As the projection is performed with the modified LuGre model, it can represent all its features and is adapted to systems with varying normal forces and arbitrary orientation of the plane of contact. The projected equations arė where G M is calculated from the global sliding velocity according to Then all components of the force can be computed from the corresponding deflection: In a steady state the bristle deflection corresponds to the forces are given by and the friction force is parallel to the tangential velocity The model yields the same steady-state characteristics as presented in Sect. 2.2: Indeed, from the steady state form we have So the sum is which can be written as a remarkable product: So the last equation can be rewritten as so that As Eqs. (29) and (37) are the same, it indicates that the steady-state behavior of the projected model is the same as that of the modified version of the LuGre model (Sect. 2.2). This set of equations can easily be implemented in multibody software. Each component of the friction force is calculated, and the resultant friction force is correctly oriented even in pure stick. These equations have exactly the same features as the modified LuGre friction model (Sect. 2.2).

Application on a simple case
To demonstrate the ability of the projected equations to be used in all directions, with multiple contact points and with varying normal force, let us consider a simple system. It consists of a cube that lies on the X-Y plane. The cube is of dimension 100 mm ×100 mm×100 mm. This cube enjoys six degrees of freedom and stands on six contact points uniformly placed on a circle of radius 50 mm (Fig. 4).
In the case of a contact between flat surfaces, it is important to identify if there is a relative rotation between the surfaces. If the parts in contact do not transmit moments (unequal distribution of normal forces along the surface), then one contact point is enough. Otherwise, at least three points are needed with four giving a more regular distribution. In this simulation, at least three contact points should be used as moments are transmitted, but for illustration, it was chosen to place six contact points to demonstrate the ability of the friction model to adapt in every direction. Therefore in this configuration the contact area is Fig. 4 Cube lying on the X-Y plane and its 6 contact points uniformly placed on a circle represented as a hexagon. In the case of more complex contact geometry, the quality of the meshing influences the results. The behavior of the cube is simulated in EasyDyn [31] using the minimal coordinate approach, and the differential equations are solved by the Newmark scheme with an automatic time step management with a limited maximum step of 0.0001 s [31]. The simulation time is set to 4 seconds. At the beginning of the simulation the cube is at rest. Only the vertical force F n is applied and is equal to 10 N. During the first second of the simulation, no other external force is applied to let the cube reach its equilibrium position. Then the pulling force P increases linearly between 1 and 2 seconds from 0 up to 0.75 N. Between 2 and 3 seconds the system stays in this state. Between 3 and 4 s the force F n starts to decrease linearly to attain zero in 4 seconds.
The friction force is computed with the projected friction model, Table 2 shows the parameters used for the modified LuGre model. These parameters are chosen as in [14]. The normal force N is computed at each contact point according to the Hunt-Crossley formula (19).
A total of 4 simulations are performed, the pulling force is in the X-Y plane and set to an angle with the X axis equal to 0 • , 30 • , 60 • , and 90 • . These 4 simulations highlight that the system is capable to reproduce pure stick and is able to withstand varying normal force under different directions of load. Figures 6-9 expose the friction force on each contact point (from F 1 to F 6) for different orientations of the pulling force P . Figure 5 shows the total friction force F tot as well. The total friction force increases up to 0.75 N, which is equal to the pulling force. As the force N decreases, at 3.5 seconds, there is a transition from stick to slip. The general behavior of the system is the same as the modified LuGre model [14]. Due to the transfer of mass, depending on the direction of the pulling force, the distribution of the normal and friction forces among the six contact points varies, but the total friction force remains the same in all cases. As the system has six contact points, the same behavior is observed every 60 • . Indeed, Figs. 6 and 7 are similar, as well as Figs. 8 and 9. This simple experiment highlights the ability of the projected equations to represent pure stick and to automatically orient the friction force in the proper direction while the system undergos varying normal force. This simulation can also be done with four contact points placed on the edges of the contact surface of the cube to represent a square contact surface. In this case, similar results are obtained, and the total friction force F tot is the same.

General kinetic friction model
The general kinetic friction (GKF) static model is a commonly used model in multibody dynamics. This model is more advanced than Coulomb's one. The GKF model has the following features: static friction, dynamic friction, Stribeck zone, and can represent the viscosity effect. The GKF model is velocity dependent, and thus pure stick cannot be represented.
To be computed numerically in multibody software, the GKF model is regularized. There are two equations: the former in "stick" regime and the latter "slip" regime [30]. Figure 10 shows the evolution of the friction force. Three zones are clearly identifiable: 1. The regularization zone from 0 to the limit velocity threshold . This zone ensures a smooth transition from −F s to F s and avoids the indetermination of the friction force at V = 0. The velocity threshold is chosen so that for the considered application, stiction is assumed when |V | < . When V is larger than , the friction force follows the evolution given by the "slip" regime of Eq. (38). 2. Between and V st , there is the Stribeck zone. 3. When V is greater than V st , the friction is in slip regime. The slope is given by the viscous coefficient f v .

Application on an experimental setup
We now apply the generalized friction model (Sect. 2.3) is now applied on a more complex multibody system problem. This system, presented in [32], represents turbines blades linked by friction elements. During the operation of a turbine, the blades are subjected to excitation that can correspond to their eigenfrequencies. To avoid fatigue and mechanical alteration of the blades, damping is necessary. This effect can be achieved by the introduction of a friction element between the blades. Out of resonance, friction is sufficient to prevent any relative motion, but close to resonance, a relative motion takes place such that the friction element acts as a damper due to the dissipation of energy by friction. Different friction elements exist; [29,33,34] study the introduction of an intermediate body between the blades. Hereafter, the friction element is directly linked to the blades in a so-called "tie-boss" coupling [32,35].
The studied system is the experimental setup investigated in [32] (Fig. 11) and consists of three identical beams coupled by a tie-boss. The excitation force is applied at the tip of the central blade (Fig. 11). Friction takes place at the contact surfaces, and these surfaces are prestressed with a normal force up to 5 N. It was observed in [32] that friction contributes to decrease the vibration amplitude, whereas all the blades preserve their single eigenmode vibration.    In the experimental setup studied in [32], the first eigenfrequency of one blade is around 50 Hz. In EasyDyn [31], similar results are obtained for one blade: the first eigenfrequency is at 53 Hz, the second at 321 Hz, and the third at 408 Hz (Table 3). Similar results were obtained by computing the eigenfrequencies of one blade under Solidworks2020 (Table 4). The software's automatic mesh generator was used and set to "high quality mesh", which generated 52878 parabolic tetrahedral solids elements and 80949 nodes. The modes were computed with the FFEPlus solver, and the first one is presented in Fig. 12. A vibration analysis of the three blades linked together through rigid connections was also realized in Solidworks2020 following the same approach, and the results are presented in Table 5 and Fig. 13. It is observed that if the blades are coupled together, then the first eigenfrequency of the system shifts to 130 Hz.
The multibody representation of a blade is a composition of several rigid and flexible bodies (Fig. 14). In Fig. 14, Body [0], Body [2] and Body [4] represent the ground, the tie-boss, and the tip of the blade, respectively. Their geometrical and inertia properties are presented in Table 6. The flexible bodies consist of two successive Euler-Bernoulli beam elements [36] with rectangular cross-section whose characteristics are presented in Table 7. The rigid bodies are linked to the corresponding nodes of the flexible bodies. Each node in the flexible bodies introduces six degrees of freedom, a blade with two elements per flexible body thus involves 24 degrees of freedom. As the system comprises three blades, the whole system involves 72 degrees of freedom. To represent the contact between two blades, one or four contact points are uniformly placed on one contact surface and interact with a plane placed on the opposite contact surface (Fig. 15). The normal force is computed at each contact point with the Hunt-Crossley formula [28]. As three differential equations are associated with each contact point (see Sect. 2.3), and there are two contact areas, there is a total of 6 or 24 supplementary differential equations to compute the friction with the projected LuGre model.
The system is simulated either with the GKF and the LuGre models. The friction parameters for the LuGre model are presented in Table 8. The parameter σ M 0 is chosen to present a sufficient stiffness (2 · 10 5 N/m with a normal force of 5 N) with respect to the stiffness of the blade in the same direction (about 2 · 10 4 N/m). The parameter σ M 1 is chosen in accordance with [10] and [14], whereas σ 2 is kept null. The limit velocity of the GKF model is equal to 0.0003 m/s, sufficiently small with respect to V st equal to 0.003 m/s.    All further simulations of this system take place in the EasyDyn framework [31]. For better control of high-frequency modes of the blades, the differential equations are integrated with the generalized α method, with a spectral radius at infinity equal to zero and a maximum time step of 10 −5 s. In all simulations, blades left and right are subjected to a crushing force perpendicularly to the contact surface to induce the contact prestress. As illustrated in Fig. 16, the latter are progressively increased during the first second of simulation up to the reach of a normal force equal to 5 N as in [32].

Simulation 1 -transition from stick to slip at low frequency
To illustrate the transition from stick to slip with both friction models, a harmonic force at 30 Hz is applied on the central blade, with a magnitude increasing linearly from 0 to 5 N during the interval [0.2-1.8] s (Fig. 17). The chosen frequency is below the first eigenfrequency of the blade alone and of the three-blade bundle.  • The transition from stick to slip takes place around instant 1 s, which corresponds to an excitation force equal to 2.5 N. The transition is smoother for LuGre than for GKF, which is explained by a rather low value of the limit velocity: indeed, in the GKF model, as far as the friction force is below the friction limit, the slip velocity is maintained below the limit velocity (which corresponds to a maximum displacement equal to 1.5 µm at 30 Hz). On the contrary, the LuGre model behaves essentially like a spring-damper system and offers more compliance with the chosen parameters. Let us mention that the order of magnitude

Simulation 2 -swept sine in the 30-80 Hz range
This second simulation concerns the frequency response of the system. In [32] the magnitude of the displacement is monitored, and the system measures the excitation force necessary to impose a displacement magnitude of 20 µm in the 30-80 Hz frequency range. In this paper, for simplicity, the force is maintained at 2 N, i.e., in the stick regime at low frequency, whereas a logarithmic swept sine is applied in the same frequency range. Figures 24 and 25 show the time history of the tip displacements during the swept sine excitation, whereas Figs. 26 and 27 represent the magnitude of the tangential velocity. In both cases, resonance takes place at about 6 s, corresponding to a frequency of about 55 Hz close to the eigenfrequency of the blade alone. It turns out that the system plays its role: blades vibrate independently of each other, and the relative slip dissipates energy by friction. At resonance, as the slip regime is well established, both models provide a similar response.
The tangential velocity exhibited by both models is also in good agreement. However, as observed in the previous simulation, the tip displacement and the tangential velocity before resonance are larger with the LuGre model. Figures 28 and 29, present the X component of   Tangential force the tangential force in terms of the X component of the tangential velocity, around time 2 s, when the excitation frequency is close to 30 Hz. In the GKF model the friction force remains in the regularization zone, and the contact behaves as a pure damper. On the contrary, in the LuGre model a hysteresis appears due to the combination of elasticity and damping. After resonance, both models exhibit comparable slip velocities obviously out of stick as larger than V st .
The evolution of the magnitude of the tangential force is plotted in Figs. 30 and 31. Forces are comparable out of resonance but not at resonance. In the GKF model, as the tangential velocity increases at resonance, the friction force decreases, accordingly with the Stribeck phenomenon. For LuGre, the dynamic nature of the model makes that, as illustrated in Fig. 32, the friction coefficient can locally take values over μ s . As a consequence, the friction force increases. Finally, Figs. 33 and 34 display the time history of the total power lost in contacts and its integral over time, i.e., the total energy loss. Again, both models give similar results. The total energy loss for LuGre (0.15 J) is a bit larger than in the GKF case (0.12 J), the difference coming essentially from the larger slips observed with LuGre before resonance.

Simulation 2 -sensitivity of the model to friction parameters
The sensitivity of the results to the contact parameters must be underlined before concluding. Concerning the GKF model, Figs. 35 and 36 show the results of the previous simulations for different values of the Stribeck velocity V st and the threshold velocity . It can be observed that changing from 0.0003 to 0.001 only has a weak impact on the results, the threshold  In the same way the influence of the damping coefficient σ 1 is provided in Figs. 37 (σ M 1 = 30 s/m) and 38 (σ M 1 = 120 s/m). Whereas a decrease from 60 to 30 does not significantly change the response of the system, the increase from 60 to 120 considerably decreases the amplification at resonance. As the purpose of the simulations is to properly tune the crushing force to get slip at the blade resonance frequency and stick otherwise, it is necessary to dispose of consolidated values of the friction parameters.

Conclusion
In this paper, we propose and apply a projected version of the modified LuGre model on a simple system and a real case problem. The projection of the modified LuGre model allows us to retrieve the direction of the friction force and is usable in a system subjected to varying normal forces. The equivalence of the projected model with the modified LuGre model is demonstrated analytically in Sect. 2.3. The proposed set of equations has all the friction features of the modified LuGre model: stick-slip, presliding displacement, frictional lag, varying break-away force, viscous friction, the Stribeck effect, and is adapted to varying normal forces. It is also adapted to be easily implemented in multibody dynamics software. This set of equations was applied on a simple system that consists of a cube that slides on a plane. The pulling force on the cube was oriented in different directions to test the friction phenomena. This simple system demonstrates the ability of the projected equations to be used in all directions and with varying normal forces. The friction model was implemented in the EasyDyn framework [31] with success and was able to obtain the same behavior as the modified LuGre model presented in [14].
To prove that this friction model can also be used in real-world applications, a system consisting of three turbine blades is presented in Sect. 5. This system investigates the use of friction to damp the response around resonance. The projected equations are used to model friction, and the results show a correlation between the experimental results presented in [32] and the results presented in this paper. In particular, the model allows us to reproduce microslip and sliding regimes depending on the excitation frequency and magnitude.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.