Paulo Flores

Margarida Machado

Miguel T Silva

Jorge M Martins

On the continuous contact force models for soft materials in multibody dynamics

Keywords: Contact force, Continuous analysis, Soft materials, Coefficient of restitution, Elastic energy, Internal damping, Multibody dynamics

A general and comprehensive analysis on the continuous contact force models for soft materials in multibody dynamics is presented throughout this work. The force models are developed based on the foundation of the Hertz law together with a hysteresis damping parameter that accounts for the energy dissipation during the contact process. In a simple way, these contact force models are based on the analysis and development of three main issues: (i) the dissipated energy associated with the coefficient of restitution that includes the balance of kinetic energy and the conservation of the linear momentum between the initial and final instant of contact; (ii) the stored elastic energy, representing part of initial kinetic energy, which is evaluated as the work done by the contact force developed during the contact process; (iii) the dissipated energy due to internal damping, which is evaluated by modeling the contact process as a single degree-of-freedom system to obtain a hysteresis damping factor. This factor takes into account the geometrical and material properties, as well as the kinematic characteristics of the contacting bodies. This approach has the great merit that can be used for contact problems involving materials with low or moderate values of coefficient of restitution and, therefore, accommodate high amount of energy dissipation. In addition, the resulting contact force model is suitable to be included into the equations of motion of a multibody system and contributes to their stable numerical resolution. A demon-strative example of application is used to provide the results that support the analysis and discussion of procedures and methodologies described in this work.

Introduction

Over the last century, the academic and scientific communities have demonstrated an increasing interest in solving problems related to contact-impact phenomena in mechanical multibody systems. Undeniably, the classical problem of the contact mechanics is a quite old topic in engineering applications [START_REF] Stolarsky | Tribology in Machine Design[END_REF][START_REF] Flores | Kinematics and Dynamics of Multibody Systems with Imperfect Joints: Models and Case Studies[END_REF]. The pioneering work on the frictionless collision between rigid bodies was due to Hertz who developed the theory with his name during the Christmas vacation of 1880 at the age of twenty-three [START_REF] Johnson | One hundred years of Hertz contact[END_REF]. The contact problem has been an intensive topic of research over the last decades and deserved the attention of many authors, which led to the development of relevant work and even to the publication of a good number of textbooks totally devoted to this issue, such as those by Goldsmith [START_REF] Goldsmith | Impact, The Theory and Physical Behaviour of Colliding Solids[END_REF], Brach [START_REF] Brach | Mechanical Impact Dynamics, Rigid Body Collisions[END_REF], Pfeiffer and Glocker [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF], Johnson [START_REF] Johnson | Contact Mechanics[END_REF], Stronge [START_REF] Stronge | Impact Mechanics[END_REF], Wriggers [START_REF] Wriggers | Computational Contact Mechanics[END_REF], among others. The contactimpact study in multibody dynamics has received a great deal of attention in the past few years and it still remains an active field of research and development [START_REF] Hippmann | An algorithm for compliant contact between complexly shaped bodies[END_REF][START_REF] Gonthier | A regularized contact model with asymmetric damping and dwell-time dependent friction[END_REF][START_REF] Flores | Dynamic analysis for planar multibody mechanical systems with lubricated joints[END_REF][START_REF] Sharf | A contact force solution for non-colliding contact dynamics simulation[END_REF][START_REF] Sousa | Development of generic multibody road vehicle models for crashworthiness[END_REF][START_REF] Djerassi | Collision with friction; Part A: Newton's hypothesis[END_REF][START_REF] Djerassi | Collision with friction; Part B: Poisson's and Stronge's hypotheses[END_REF][START_REF] Bowling | Energetically consistent simulation of simultaneous impacts and contacts in multibody systems with friction[END_REF][START_REF] Dopico | Dealing with multiple contacts in a human-in-theloop application[END_REF][START_REF] Ambrósio | Improved bushing models for general multibody systems and vehicle dynamics[END_REF][START_REF] Mukras | Analysis of planar multibody systems with revolute joint wear[END_REF][START_REF] Choi | An efficient and robust contact algorithm for a compliant contact force model between bodies of complex geometry[END_REF][START_REF] Flores | Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach[END_REF].

Contact-impact events can frequently occur in multibody systems, and in many cases the function of mechanical systems is based on them. Gilardi and Sharf [START_REF] Gilardi | Literature survey of contact dynamics modeling[END_REF] describe an impact as a complex physical phenomenon for which the main characteristics are a very short duration, high force levels, rapid energy dissipation, and large changes in the velocities of bodies. In turn, contact implies a continuous process which takes place over a finite time. The analysis of contact between two bodies can be extended to the analysis of collision in a multibody system. Whenever a body in a multibody system experiences a hard stop, contact forces of complex nature act on the body and the corresponding impulse is transmitted throughout the system.

The bodies that belong to a multibody system can be considered as rigid or flexible. A body is said to be rigid when its deformations are assumed to be small such that they do not influence the global motion produced by the body. The expression flexible multibody system refers to a system holding deformable bodies with internal dynamics [START_REF] Shabana | Dynamics of Multibody Systems[END_REF]. In fact, rigid bodies are a representation of reality because bodies are not absolutely rigid in nature. However, many physical mechanical systems contain a combination of rigid and flexible bodies, or rigid bodies with soft surfaces, or bodies that are rigid enough to be considered rigid overall, but soft enough to experience significant local deformations during contact. In this case, the contact approach must be able to model the dynamics of contact between compliant surfaces.

The compliant continuous contact force models, commonly referred as penalty methods, have been gaining significant importance in the context of multibody systems with contacts due to their computational simplicity and efficiency. In these models, the contact force is expressed as a continuous function of the relative deformation of the contacting surfaces during the complete period of contact. In addition, most of the continuous contact force models contribute to an efficient integration of the equations of motion and account for some level of energy dissipation. This amount of good reasons explains why some of the most popular commercial codes use the penalty methods to treat contact problems [START_REF] Ryan | ADAMS-Multibody System Analysis Software, Multibody Systems Handbook[END_REF][START_REF] Smith | DADS-Dynamic Analysis and Design System, Multibody Systems Handbook[END_REF][START_REF]Visual NASTRAN 4D[END_REF].

The energy loss due to the contact process in the normal direction can be expressed in terms of the coefficient of restitution. This parameter has different definitions, being one of the most popular and commonly used the Newton's law of restitution, also known as kinematical coefficient of restitution. In a simple manner, this definition states that when two bodies collide, the post-impact velocity is related to the preimpact velocity by a constant of proportionality called coefficient of restitution [START_REF] Lee | On the dynamics of intermittent-motion mechanisms, part 1: dynamic model and response[END_REF]. This coefficient is, in general, assumed to be constant; however, it is dependent on several factors such as the geometry of the contacting surfaces, preimpact velocity, local material properties, duration of the contact, temperature, and friction [START_REF] Goldsmith | Impact, The Theory and Physical Behaviour of Colliding Solids[END_REF][START_REF] Brach | Mechanical Impact Dynamics, Rigid Body Collisions[END_REF][START_REF] Jackson | Predicting the coefficient of restitution of impacting elasticperfectly plastic spheres[END_REF]. The analysis and discussion of these issues is out of the scope of the present work, being the interested reader in this particular topic referred to the work by Gilardi and Sharf [START_REF] Gilardi | Literature survey of contact dynamics modeling[END_REF].

The initial kinetic energy produced in a contact event is converted in different measures during the contact process. First of all, part of the initial kinetic energy is transformed into elastic strain energy which, in turn, is converted in the final kinetic energy at the end of the contact. The remaining amount of kinetic energy is dissipated during the contact process in several forms such as wave propagation [START_REF] Stronge | Impact Mechanics[END_REF][START_REF] Barkan | Impact design. In: Mechanical Design and Systems Handbook[END_REF], material damping [START_REF] Goldsmith | Impact, The Theory and Physical Behaviour of Colliding Solids[END_REF][START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF], plasticity [START_REF] Brach | Mechanical Impact Dynamics, Rigid Body Collisions[END_REF][START_REF] Lankarani | Continuous contact force models for impact analysis in multibody systems[END_REF], and heat [START_REF] Gilardi | Literature survey of contact dynamics modeling[END_REF][START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF]. Stronge [START_REF] Stronge | Impact Mechanics[END_REF] demonstrated that stress and deformation waves do not contribute to energy dissipation when the size and material properties of the contacting bodies are similar. Shivaswamy [START_REF] Shivaswamy | Modeling contact forces and energy dissipation during impact in multibody mechanical systems[END_REF] showed that at low impact velocities, the damping force component is the prime factor for energy dissipation. At high impact velocities, that is, those exceeding the propagation velocity of elastic deformation waves are likely to dissipate energy in a form of permanent local deformation.

The general motivation for this work comes from current interest in developing mathematical models for the dynamics of constrained multibody systems involving contact events between soft materials, such as in the case of biomechanical systems and bushing elements made with polymer materials [START_REF] Machado | Development of a planar multibody model of the human knee joint[END_REF][START_REF] Meireles | Strain shielding in distal femur after patellofemoral arthroplasty under different activity conditions[END_REF][START_REF] Bei | Multibody dynamic simulation of knee contact mechanics[END_REF][START_REF] Lin | Simultaneous prediction of muscle and contact forces in the knee during gait[END_REF][START_REF] Burgin | Impact testing to determine the mechanical properties of articular cartilage in isolation and on bone[END_REF][START_REF] Piazza | Three-dimensional dynamic simulation of total knee replacement motion during a step-up task[END_REF][START_REF] Silva | Evaluation of the contact forces developed in the lower limb/orthosis interface for comfort design[END_REF]. In fact, modeling contacts with a high degree of energy dissipation plays a crucial role in the analysis, design and control of many mechanical systems in the measure that soft materials are characterized by low or medium values of the coefficient of restitution. The methodology adopted in this work follows closely that of [START_REF] Lankarani | Canonical equations of motion and estimation of parameters in the analysis of impact problems[END_REF][START_REF] Ye | A note on the Hertz contact model with nonlinear damping for pounding simulation[END_REF], in which an explicit relation between the coefficient of restitution and a hysteresis damping factor is derived. Furthermore, the approach presented in this study is based on the analysis of dissipated energy related to the coefficient of restitution, stored elastic energy, and the dissipated energy associated with the internal damping of the contacting bodies. The methodology and procedures adopted throughout this work are presented and analyzed with the help of a numerical simulation of the classical bouncing ball problem.

General issues of continuous contact force models

The theory developed by Hertz remains the foundation for most contact problems encountered in engineering applications. The Hertz law relates the contact force with a nonlinear power function of deformation and is expressed as [START_REF] Hertz | On the contact of solids-on the contact of rigid elastic solids and on hardness[END_REF]

F N = Kδ n (1)
where K represents the generalized stiffness parameter and δ is the relative normal deformation between the contacting bodies. The exponent n is typically equal to 3/2 [START_REF] Flores | Kinematics and Dynamics of Multibody Systems with Imperfect Joints: Models and Case Studies[END_REF]. The generalized parameter K is dependent on the material properties and the shape of the surfaces in contact. For two spheres i and j in contact, the generalized parameter is a function of radii of the spheres and the material properties as [START_REF] Goldsmith | Impact, The Theory and Physical Behaviour of Colliding Solids[END_REF]

K = 4 3(σ i + σ j ) R i R j R i + R j 1 2 (2)
in which the material parameters σ i and σ j are given by

σ l = 1 -υ 2 l E l , (l = i, j ) (3)
and the quantities υ l and E l are the Poisson's ratio and Young's modulus associated with each sphere, respectively. It is apparent that the Hertz contact law given by ( 1) is limited to contact with elastic deformations and does not include the damping effect. This force model represents the contact process as a nonlinear spring along the direction of collision. Although the Hertz law is based on the elasticity theory [45], some studies have been performed to extend its application to include energy dissipation. An alternative contact force approach that accounts for the energy loss during the contact process is the Kelvin-Voigt model. This force model considers a linear spring in conjunction with a linear damper and can be expressed as [START_REF] Khulief | A continuous force model for the impact analysis of flexible multibody systems[END_REF]

F N = Kδ + D δ (4)
where the first term of the right-hand side is referred to as the elastic force and the second term accounts for the energy dissipation during the contact process. In (4), D is a hysteresis coefficient and δ represents the relative normal contact velocity, having the remaining variables the same meaning as it was described above. The constitutive law given by (4) presents some weaknesses, namely the fact that the contact force at zero deformation is not continuous due to the existence of the damping component. This particular issue is not realistic because when the contact begins, both elastic and damping force components must be null. Moreover, at the end of the contact, the deformation is null, the relative contact velocity is negative and, consequently, the resulting contact force is also negative, which does not make sense from the physical point of view because the bodies can not attract each other [START_REF] Johnson | One hundred years of Hertz contact[END_REF][START_REF] Marhefka | A compliant contact model with nonlinear damping for simulation of robotic systems[END_REF][START_REF] Bibalan | A study of soft contact models in simulink[END_REF].

Hunt and Crossley [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF] showed that the linear spring-damper model does not represent the physical nature of energy transferred during the contact. Instead, they represented the contact force by the Hertz force-deformation law with a nonlinear viscous-elastic element as

F N = Kδ n + χδ n δ ( 5
)
where χ is called the hysteresis damping factor given by

χ = 3K(1 -c r ) 2 δ(-) (6)
in which K represents the generalized stiffness parameter, c r is the coefficient of restitution and δ(-) is the initial contact velocity. After some mathematical manipulation, the expression for this contact force model has the following form [START_REF] Dopico | Dealing with multiple contacts in a human-in-theloop application[END_REF]:

F N = Kδ n 1 + 3(1 -c r ) 2 δ δ(-) (7)
With this model, the energy loss during the contact is assumed to be associated with the material damping of the contacting bodies, which would dissipate energy in the form of heat. The Hunt and Crossley force model expresses the damping as a function of deformation, which sounds physically reasonable. Furthermore, this model does not present discontinuities at the initial instant of contact and at the end of contact, that is, it starts and ends with the zero value.

Lankarani and Nikravesh [START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF] used the general trend of the Hertz contact law incorporated with a hysteresis damping factor to propose a continuous force model for the contact analysis of multibody systems. The hysteresis damping factor that accommodates the energy dissipation during the contact is a function of the deformation and is written as

χ = 3K(1 -c 2 r ) 4 δ(-) (8)
Based on the analysis of Hunt and Crossley [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF] that separates the contact force into elastic and dissipative components, the force model proposed by Lankarani and Nikravesh can be expressed as

F N = Kδ n 1 + 3(1 -c 2 r ) 4 δ δ(-) (9)
which is quite similar to the model proposed by Hunt and Crossley. It should be highlighted that in the Lankarani and Nikravesh force model, the parameters have been evaluated from the physics of the problem by equating the kinetic energy loss in the system to the work done by the contact force. Since the kinetic energy loss is in terms of the relative velocity squared, then the expression for the hysteresis damping factor and the contact force include a term in c 2 r . The contact force models given by ( 7) and ( 9) are valid for the cases in which the dissipated energy during the contact is small when compared to the maximum absorbed elastic energy, that is, the relation is valid for the values of the coefficient of restitution close to unity. The consequences of this simplifying premise are discussed in detail by Lankarani and Nikravesh [START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF]. In a later paper, Lankarani and Nikravesh [START_REF] Lankarani | Continuous contact force models for impact analysis in multibody systems[END_REF] proposed a new approach for the contact force model, in which the permanent or plastic deformation was also included.

By analyzing ( 6) and ( 8), it can be observed that for a perfectly elastic contact, i.e., c r = 1, the hysteresis damping factor assumes a zero value, while for a perfectly plastic contact, i.e., c r = 0, the hysteresis damping factor does not assume an infinite value as it was expected. This fact is not surprising because, for instance, Lankarani and Nikravesh [START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF] derived their model for high values of restitution coefficient, being, therefore, valid for hard contacts, such as those between metals [START_REF] Flores | Modeling and simulation of wear in revolute clearance joints in multibody systems[END_REF][START_REF] Tian | Dynamics of spatial flexible multibody systems with clearance and lubricated spherical joints[END_REF][START_REF] Flores | Influence of the contact-impact force model on the dynamic response of multibody systems[END_REF][START_REF] Tian | Simulation of planar flexible multibody systems with clearance and lubricated revolute joints[END_REF]. The values of the coefficient of restitution for soft materials are low or medium [START_REF] Piazza | Three-dimensional dynamic simulation of total knee replacement motion during a step-up task[END_REF][START_REF] Ye | A note on the Hertz contact model with nonlinear damping for pounding simulation[END_REF][START_REF] Bibalan | A study of soft contact models in simulink[END_REF][START_REF] Ambrósio | Sensitivity of a vehicle ride to the suspension bushing characteristics[END_REF], which clearly indicate that the hysteresis damping factor given by ( 6) and ( 8) cannot be used. Therefore, a more accurate expression for this factor is required.

Description of the elastic contact between two solid spheres

Prior to present the detailed description of the presented approach for the normal contact force, let us first introduce some fundamental issues related to contact mechanics, which play a crucial role since they constitute the basis that supports the general methodology used to develop the main procedures followed in the present work.

Among the different types of contact that can be defined for the collision between two bodies that belong to a multibody system, the one-dimensional direct central contact between two solid spheres is considered here with the purpose to keep the analysis simple.

When the spheres contact each other, deformation takes place in the local contact area resulting in a contact force. Two different phases can be distinguished during the contact process, namely the compression, approaching or loading period, and the restitution, separating or unloading period [START_REF] Goldsmith | Impact, The Theory and Physical Behaviour of Colliding Solids[END_REF][START_REF] Brach | Mechanical Impact Dynamics, Rigid Body Collisions[END_REF]. The first phase starts when the two spheres come in contact and ends when the maximum deformation is reached. During this phase, the relative velocity of the contact points on the two bodies in the normal contact direction is gradually reduced to zero. The end of the compression phase is referred to as the instant of maximum compression or maximum approach, which corresponds to the maximum normal contact force if no dissipation of energy is considered. In turn, the restitution phase starts at this instant of time and ends when the two spheres separate from each other.

Figure 1 illustrates the two different phases of a direct central contact between two solid spheres. The spheres are labeled as bodies i and j and have masses m i and m j . The pre-and post-impact velocities associated with each sphere are denoted as v (-) i , v (-) j , v (+) i , and v (+) j , respectively. At the end of the compression phase, the two spheres move with a common instantaneous velocity represented by v (max) ij .

Figure 2 shows a general diagram for the deformation evolution with time during a single contact between two spheres, as well as the contact force plotted against deformation. In these diagrams, the instant of initial contact, the instant of maximum deformation and the instant of separation are denoted as t (-) , t (max) , and t (+) , respectively. In turn, the symbols δ max and F max refer to the values of deformation and contact force at the instant of time t (max) , respectively.

Dissipated energy associated with the coefficient of restitution

The simplest manner to quantify the energy loss during a contact event is to use the concept of coefficient of restitution, which can be evaluated from the balance of energy that occurs between the initial instant of contact t (-) and the final instant of contact t (+) . For this purpose, let consider T (-) and T (+) , respectively, the kinetic energy of the two contacting spheres at the instants of time t (-) and t (+) . Thus, the energy balance can be expressed as [START_REF] Beer | Vector mechanics for engineers[END_REF]

E = T (-) -T (+) = 1 2 m i v (-) i 2 -v (+) i 2 + 1 2 m j v (-) j 2 -v (+) j 2 ( 10
)
in which, the post-impact velocities are the unknown quantities.

The linear momentum associated with the two spheres i and j in the direction of the contact is conserved because the contact forces are internal to the system [START_REF] Greenwood | Principles of Dynamics[END_REF]. Therefore, it is possible to write the following equation:

m i v (-) i -v (+) i + m j v (-) j -v (+) j = 0 (11)
On the other hand, taking advantage of the kinematical or Newton's definition of the coefficient of restitution, which expresses the ratio between the post-and the pre-impact velocities of the spheres in the contact, yields [START_REF] Stronge | Impact Mechanics[END_REF]

c r = - v (+) i -v (+) j v (-) i -v (-) j = - δ(+) δ(-) (12)
where δ(-) and δ(+) are, respectively, the relative approaching velocity and the relative separating velocity. Equations ( 10)-(12) form a system of three equations with three unknowns v (+) i , v (+) j , and E. Solving this system of equations results in

v (+) i = 1 m i + m j (m i -c r m j )v (-) i + (1 + c r )m j v (-) j ( 13
)
v (+) j = 1 m i + m j (m j -c r m i )v (-) j + (1 + c r )m i v (-) i (14) E = 1 2 m i m j m i + m j v (-) i -v (-) j 2 1 -c 2 r ( 15
)
Equation ( 15) represents the dissipated energy in the contact process of two solid spheres in terms of the coefficient of restitution and relative approaching velocity, which can be simplified in the following form:

E = 1 2 m δ(-) 2 1 -c 2 r ( 16
)
where m is the equivalent mass of the two spheres given by m = m i m j m i + m j [START_REF] Bowling | Energetically consistent simulation of simultaneous impacts and contacts in multibody systems with friction[END_REF] The coefficient of restitution describes the type of collision. For a fully elastic contact the restitution coefficient is equal to unity, while for a fully plastic contact the restitution coefficient is null. The most general and predominant type of collision is the oblique eccentric collision, which involves both relative normal velocity and relative tangential velocity [START_REF] Maw | The oblique impact of elastic spheres[END_REF][START_REF] Zukas | Impact Dynamics[END_REF].

Stored elastic strain energy during contact

As it was mentioned above, at the end of the compression phase, the two spheres move with a common instantaneous velocity v (max) ij . Moreover, the initial kinetic energy associated with the contact of the two spheres has three different destinations. One part is transformed into the kinetic energy of the two spheres moving with the same common velocity. A second part is converted to stored elastic strain energy. Finally, a third part of the initial kinetic energy is dissipated. Thus, the balance of energy between the start and end of the compression phase can be written as follows:

T (-) = T (max) + U (max) + E c ( 18
)
in which T (max) represents the kinetic energy of the system at the end of the compression phase, U (max) is the maximum elastic strain energy stored and E c denotes the dissipated energy associated with internal damping of the material. This last term will be discussed in detail in the next section.

The stored strain energy is equal to the work done by the contact force that develops from the state of zero deformation to the state of maximum deformation, that is [START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF]

U (max) = δ max 0 F N dδ (19)
In order to obtain a solution of this integral let consider that the contact force is represented by the pure elastic Hertz law (1), and consequently, the elastic strain energy absorbed can be rewritten as

U (max) = δ max 0 Kδ 3 2 dδ (20)
Finally, the stored elastic strain energy is given by

U (max) = 2 5 Kδ 5 2 max ( 21
)

Dissipated energy due to internal damping

In this section, the dissipated energy due to internal damping is determined. With the purpose of reaching this desideratum, let consider an equivalent system to the contact process between two spheres materialized by a single degree of freedom dynamic system, as depicted in Fig. 3. In this model, the initial deformation between the spheres is zero and the relative initial contact velocity is denoted by δ(-) . The mathematical representation of this dynamic system can be expressed as [START_REF] Brach | Mechanical Impact Dynamics, Rigid Body Collisions[END_REF][START_REF] Hartog | Mechanical Vibrations[END_REF][START_REF] Steidel | An Introduction to Mechanical Vibrations, 3rd edn[END_REF]]

m δ + c δ + kδ = 0 ( 22
)
where m is the equivalent mass given by ( 17), c represents the viscous damping coefficient and k is the equivalent stiffness. The main reason to use the model of ( 22) owes to the fact that it matches precisely with the ideal Kelvin-Voigt viscoelastic approach [START_REF] Goldsmith | Impact, The Theory and Physical Behaviour of Colliding Solids[END_REF]. Equation ( 22) is a second-order differential equation with constant coefficients; the solution of which is given by [START_REF] Hartog | Mechanical Vibrations[END_REF] δ = e -ζ ωt δ(-)

ω d sin(ω d t) (23) δ = e -ζ ωt δ(-) cos(ω d t) -ζ e -ζ ωt ω δ(-) ω d sin(ω d t) ( 24
)
where ζ denotes the damping ratio, ω is the radial frequency and ω d represents the damped frequency. These three parameters are expressed as

ζ = c 2mω (25) ω = k m ( 26
)
ω d = ω 1 -ζ 2 (27)
In the case that the damping effect is neglected, the values of t (max) and δ max are given by

δ max = δ(-) ω ( 28
)
t (max) = π 2ω (29)
Substituting ( 28) and ( 29) into ( 23) and ( 24), and taking advantage of the fundamental law of the trigonometry, yields the following mathematical relation:

δ δ max 2 + δ δ(-) 2 = 1 ( 30
)
which is an equation of an ellipse for the deformation and velocity. Thus, it is possible to express the velocity of deformation as function of the deformation for both compression and restitution phases, respectively, as [START_REF] Lankarani | Continuous contact force models for impact analysis in multibody systems[END_REF] At this stage, it is important to note that the post-impact velocity δ(+) is negative because the spheres are separated from each other.

δ = δ(-) 1 - δ δ max 2 (31) δ = δ(+) 1 - δ δ max 2
The dissipated energy is evaluated through the work done by the damping force component, which can be written as [START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF]

E = χδ 3 2 δ dδ ( 33
)
where refers to the integral around a hysteresis loop for the contact force [START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF]. The total dissipated energy due to the damping force associated with the compression and restitution phases can be expressed, respectively, as

E c = δ max 0 χδ 3 2 δ(-) 1 - δ δ max 2 dδ ( 34
) E r = δ max 0 χδ 3 2 δ(+) 1 - δ δ max 2 dδ ( 35
)
and consequently, the total amount of dissipated energy is given by

E = E c + E r = χ δ(-) + δ(+) δ max 0 δ 3 2 1 - δ δ max 2 dδ ( 36
)
Applying the integral transformation to (36) yields [START_REF] Ye | A note on the Hertz contact model with nonlinear damping for pounding simulation[END_REF]

E = χ δ(-) + δ(+) δ 5 2 max 1 0 x 3 2 1 -x 2 dx ( 37
)
It should be highlighted that the value of the integral term in [START_REF] Bei | Multibody dynamic simulation of knee contact mechanics[END_REF] can only be evaluated numerically, yielding

1 0 x 3 2 1 -x 2 dx 1 4 ( 38
)
The detailed demonstration of this solution is offered in the appendix presented at the end of this paper. Hence, substituting this value into [START_REF] Bei | Multibody dynamic simulation of knee contact mechanics[END_REF] results in

E = 1 4 χ δ(-) + δ(+) δ 5 2 max ( 39
)
Regarding now [START_REF] Burgin | Impact testing to determine the mechanical properties of articular cartilage in isolation and on bone[END_REF] and taking into account the definition of the coefficient of restitution [START_REF] Flores | Dynamic analysis for planar multibody mechanical systems with lubricated joints[END_REF], then (39) can be rewritten in the following form:

E = 1 4 χ(1 -c r ) δ(-) δ 5 2 max ( 40
)
When the dissipated energy during the contact of two spheres is taken into account, the balance of energy given by ( 18) can be expanded as

1 2 m i v (-) i 2 + 1 2 m j v (-) j 2 = 1 2 (m i + m j )v (max) ij + 2 5 Kδ 5 2 max + 1 4 χ δ(-) δ 5 2 max (41)
The linear momentum balance between instants t (-) and t (max) is given by

m i v (-) i + m j v (-) j = (m i + m j )v (max) ij ( 42
)
Thus, combining ( 16), ( 40), [START_REF] Silva | Evaluation of the contact forces developed in the lower limb/orthosis interface for comfort design[END_REF], and (42), after some basic mathematical manipulation, the hysteresis damping factor χ can finally be expressed as

χ = 8K(1 -c r ) 5c r δ(-) (43)
which represents an explicit relation of the equivalent damping factor, the generalized stiffness, the initial contact velocity, and the coefficient of restitution.

From the analysis of ( 43), it can be concluded that for a perfectly elastic contact, i.e. c r = 1, the hysteresis damping factor assumes a zero value, and when the contact is purely plastic, i.e., c r = 0, the hysteresis damping factor is infinite, which is reasonable from physical point of view. Finally, introducing the expression of the hysteresis damping factor given by ( 43) into (5) yields the expression for normal contact force described in this work

F N = Kδ n 1 + 8(1 -c r ) 5c r δ δ(-) (44)

Demonstrative application to a bouncing ball

The example of application considered here is the classical bouncing ball problem, which is one of the simplest contact systems [START_REF] Ye | A note on the Hertz contact model with nonlinear damping for pounding simulation[END_REF][START_REF] Flores | Contact-impact analysis in multibody systems based on the nonsmooth dynamics approach[END_REF]. Figure 4 shows an elastic ball made of PTFE, with an initial height equal to 1.0 m, a mass of 1.0 kg, a moment of inertia equal to 0.1 kg m 2 , a radius R = 0.1 m, and an equivalent stiffness equal to 140 × 10 6 N/m 3/2 . The ball is released from the initial position under the action of gravity only, which is taken to act in the negative y direction. Thus, the ball falls down until it collides with the ground, which is considered to be rigid and stationary. When the ball collides with the ground, a contact takes place and the ball rebounds, producing jumps, in which height depends on the value of the coefficient of restitution.

When the ball contact with the ground, the deformation is evaluated as

δ = y b -R (45)
where y b is the y coordinate of the ball center of mass and R represents the ball radius.

Positive values of δ mean that there is no contact between the ball and ground. Therefore, the detection of the instant of contact occurs when the sign of deformation changes between two discrete instants in time [START_REF] Flores | On the contact detection for contact-impact analysis in multibody systems[END_REF]. This problem was analyzed by using a FORTRAN code named MUBODYNA (acronym for Multibody Dynamics) [START_REF] Flores | MUBODYNA-A FORTRAN program for dynamic analysis of planar multibody systems[END_REF]. Figure 5 includes the global results of the ball kinematics, namely the ball position and velocity, when the coefficient of restitution is equal to 0.9 and the contact force model described in this paper is considered. When the ball hits the ground, a contact occurs and the ball is rebounded. The restitution coefficient is less than unit, therefore, the ball jumps less high after each contact. The time lapse between the contacts tends to zero and ends up when the ball remains on the ground.

Figures 6a-c show the time plots of the deformation, velocity of deformation and contact force developed during the first contact between the ball and the ground. Additionally, Fig. 6d illustrates the contact force plotted against the deformation. In the case studies presented here, five different values for the coefficient of restitution are selected, namely 1.0, 0.8, 0.6, 0.4, and 0.2. As it was expected, when the coefficient of restitution decreases, the maximum deformation and the maximum contact force values are reduced. These facts are accompanied with the increase of the contact duration and the hysteresis loop, which in turn is associated with the energy dissipated during the contact process. When the coefficient of restitution is equal to unit, which corresponds to the pure Hertz contact law, there is no energy dissipation in the contact process. This fact is evident in the force-deformation diagram of Fig. 6d, which does not present a hysteresis loop. In general, the global data produced here are corroborated by similar analysis available in the thematic literature [START_REF] Gonthier | A regularized contact model with asymmetric damping and dwell-time dependent friction[END_REF][START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF][START_REF] Ye | A note on the Hertz contact model with nonlinear damping for pounding simulation[END_REF].

A similar analysis presented by Lankarani and Nikravesh [START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF] is also included in the present study, in which the post-impact velocity that results from the first impact is used to compute the amount of restitution produced by the contact force. Thus, the post-restitution coefficient is evaluated as which differs from the coefficient of restitution provided as input to the contact force models hereafter referred as pre-restitution. Figure 7 illustrates the plots of the post and prerestitution coefficients for the three contact force models considered in the present work. The straight line that corresponds to the ratio between post-and prerestitution coefficients equal to the unity can be used as reference. From the plots of Fig. 7, it can be drawn that the Hunt and Crossley and the Lankarani and Nikravesh force models are valid for high values of the coefficient of restitution, as it was stated in reference [START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF]. Furthermore, it can also be observed that when the coefficient of restitution is null, these two force models do not reflect a perfectly plastic collision. Conversely, the contact force model described in the present work exhibits close approximation for the entire interval between 0 and 1.

c post r = - δ(+) measured δ(-) (46)
In order to better understand what are the consequences in terms of the dynamic system response when different contact force models are used, let consider a particular bouncing problem in which the ball falls down due to gravity action and impacts with a table. The table has a mass equal to 1.0 kg and is animated by a harmonic excitation given by y t (t) = 0.5 + 0.1 sin 7.5t [START_REF] Marhefka | A compliant contact model with nonlinear damping for simulation of robotic systems[END_REF] The occurrence of contact between ball and table is determined by evaluating the relative penetration at any time during the numerical solution of the dynamic equation of motion as

δ = y b -y t -R ( 48
)
where y b and y t represent the y coordinate of the ball and table centers of mass, respectively, and R is the radius of the ball. In a simple way, the ball and table approach each other with a nonzero relative velocity in the normal direction until the distance between them has been vanished. Then two different scenarios can occur, namely: (i) the ball and table separate immediately after the collision with a finite positive normal velocity, (ii) or the ball and table may remain in permanent contact. The system response scenario depends on the constitutive contact law used to model the contact event, on the value of restitution coefficient and on the overall system dynamics.

Figure 8 shows the ball motion when the contact is modeled with the pure Hertz law given by [START_REF] Stolarsky | Tribology in Machine Design[END_REF]. By analyzing the system response, it can be observed that the ball rebounds after each impact with the table. When the impact occurs in the ascendant phase of the table motion, the ball gains energy and reaches a higher position than the one before the impact.

The influence of the use of three different dissipative contact force models on the ball motion is illustrated in the plots of Fig. 9. Three difference values for the coefficient of restitution are considered, namely 0.9, 0.6, and 0.3. From the plots of Fig. 9, it can be concluded that the ball losses energy in each impact with the table. Furthermore, after some finite number of impacts the ball remains in contact with the table. This observation is valid for all dissipative contact force models used. For high values of the restitution coefficient, the Hunt and Crossley model exhibits a response similar to the case simulated with the Lankarani and Nikravesh model. The system response presents less and lower intensity rebounds when the described force approach is considered, being the permanent contact between the ball and table reached faster. This effect is more evident for lower values of the coefficient of restitution, as it can be observed in Fig. 9c. In general, it can be drawn that the Hunt and In conclusion, for a perfectly elastic impact (c r = 1), the three dissipative contact force models used reduce to the general Hertzian elastic model given by [START_REF] Stolarsky | Tribology in Machine Design[END_REF], because the hysteresis damping factor assumes a zero value. However, for a perfectly plastic impact (c r = 0), the hysteresis damping factor given by Hunt and Crossley and Lankarani and Nikravesh force models does not assume an infinite value. In sharp contrast, when the described model is considered, the hysteresis damping factor tends to infinity and, consequently, the normal contact force tends also to infinity. This case makes sense from physical point of view, in the measure that for perfectly plastic impacts, the bodies will be in permanent contact after the collision, which means that the force necessary to separate them is infinite, as predicted by the force model described throughout this work. Finally, an experimental work, such as the one presented by Shivaswamy [START_REF] Shivaswamy | Modeling contact forces and energy dissipation during impact in multibody mechanical systems[END_REF], must be developed in order to validate the model described here.

Conclusions

The problem of continuous contact force models for soft materials in multibody dynamics has been analyzed throughout this work. In the sequel of the process, the main issues associated with the fundamental contact mechanics were revised. Furthermore, some of the limitations and drawbacks of the classical contact force models commonly used when modeling contact events in multibody systems and robotics were discussed, also. In particular, the case in which the restitution coefficient exhibits low or medium value was analyzed. With the intent to overcome those drawbacks, a new approach for a continuous contact force was described, resulting in a mathematical expression that can be applied in the entire range of the coefficient of restitution, that is, in the cases of soft and hard contacts. This approach is straightforward to be implement in any multibody code for dynamic analysis, expresses a dependence on the material, mechanical, geometrical, and kinematical properties of the contacting bodies. In addition, it contributes to the stable numerical resolution of the equations of motion of constrained multibody systems. Finally, a simple problem was used to analyze and compared the different continuous contact force models includes in the present study, where the limits of application of each model was shown.

Fig. 1 Fig. 2

 12 Fig. 1 Representation of the one-dimensional direct central collision between two solid spheres: (a) Beginning of the compression phase; (b) Maximum approach instant; (c) Termination of the restitution phase

Fig. 3

 3 Fig. 3 (a) Contact between two rigid spheres; (b) Equivalent system

Fig. 4 Fig. 5

 45 Fig. 4 Bouncing ball example

Fig. 6 Fig. 7

 67 Fig. 6 Influence of the coefficient of restitution on the bouncing ball behavior: (a) Deformation versus time; (b) Velocity of deformation versus time; (c) Contact force versus time; (d) Contact force-deformation relation

Fig. 8 Fig. 9

 89 Fig. 8 Bouncing ball response for the Hertz contact force model

Acknowledgements

The authors would like to thank the Portuguese Foundation for Science and Technology (FCT) for the support given through projects PROPAFE-Design and Development of a Patello-Femoral Prosthesis (PTDC/EME-PME/67687/2006), DACHOR-Multibody Dynamics and Control of Hybrid Active Orthoses (MIT-Pt/BSHHMS/0042/2008) and BIOJOINTS-Development of advanced biological joint models for human locomotion biomechanics (PTDC/EME-PME/099764/2008). The authors want also to express their sincere gratitude to Professor Rui Pereira, from Department of Mathematics of the University of Minho, for his help in the resolution of some mathematical and numerical problems included in the present work. The second author expresses her gratitude to FCT for the Ph.D. grant SFRH/BD/40164/2007.

Appendix: Integral solution of the dissipated energy due to internal damping

This appendix contains the detailed resolution of the integral that defines the amount of energy dissipated during a contact event in terms of internal damping. For this purpose, let us rewrite (36)

Now, a substitution of variable must be done according to the following relation:

Then the lower and upper limits of integration can be obtained as follows:

Taking advantage of (A.2)-(A.5), (A.1) can be rewritten in the following form:

in which the integral term is of differential binomial type [START_REF] Piskounov | Cálculo Diferencial e Integral[END_REF] x a α + βx b c dx (A.7)

where a, b, c, α, and β are given constants. In the case of (A.6), these constants assume the following values a = 3/2, b = 2, c = 1/2, α = 1, and β = -1.

According to Piskounov [START_REF] Piskounov | Cálculo Diferencial e Integral[END_REF], the integral given by (A.7) together with the above values for a, b, c, α, and β cannot be solved analytically because it is not possible to express the integral in terms of elementary functions. Consequently, a numerical procedure must be adopted. In the present work, the trapezoidal rule is chosen to obtain the solution of this integral, which can be written as [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF] x n x 0 y(x) dx

x nx 0 2n (y 0 + 2y 1 + 2y Finally, the result of the integral that appears in (A.6) is obtained by performing the computational simulation of the above procedure, yielding area = 0.24971047130972 1/4.


Götz Von Groll

David J Ewins

The harmonic balance method with arc-length continuation in rotor/stator contact problems

There are a variety of abnormal running conditions in rotating machinery which lead to rotor/stator interaction dynamics which, in turn, can cause a rich mixture of effects associated with rub-related phenomena. These effects manifest themselves in the occurrence of multiple solutions for steady-state vibration response scenarios, including amplitude jumps during rotor acceleration, and vibration responses at different/multiple frequencies of excitation forces such as imbalance. This paper describes a numerical algorithm based on the harmonic balance method to calculate the periodic response of a non-linear system under periodic excitation. The algorithm also calculates the stability of the periodic solutions found, marks turning and bifurcation points, and follows a solution branch over varying system parameters via arc-length continuation.

INTRODUCTION

The motivation for this study comes from rotor/stator contact induced vibration problems in turbo-machinery. They can include: rotors touching seals, rotor touching retainer bearings when main active magnetic bearings fail, inter-shaft contact in multiple spool engines, rotor blades contacting the stator, increased bearing clearance through wear or outright bearing failure. In many of these scenarios the rotor continues to rotate and so, depending on the problem, it is often the steady-state response to the out-of-balance excitation forces which is of concern, rather than a particular transient event.

Both the harmonic balance method (HBM) and continuation schemes are wellknown numerical tools to study nonlinear dynamics problems. However, they seem to be used rarely in conjunction with each other in engineering applications, as continuation appears more frequently with time-domain methods, such as shooting or boundary value problem solvers. Recent rotor dynamics examples using timedomain methods with continuation can be found in Sundararajan and Noah [START_REF] Sundararajan | Transactions of the ASME[END_REF]2] dealing with squeeze-film-damper and journal bearing analysis and Petrov [START_REF] Petrov | Analysis of periodic regimes of forced[END_REF] for shroud/blade friction. Contrary to the shooting method or boundary value prob-lems, the HBM is essentially a frequency-domain method. It is not only convenient for purposes of linearisation of systems with small non-linearities [START_REF] Géradin | Fourth International Conference on Vibrations in Rotating Machinery[END_REF], but can can also be applied to large non-linearities [START_REF] Choi | [END_REF]. Kim et al. [6] analysed the behaviour of rotors with bearing housing clearances using the HBM at discrete speeds but without continuation.

As with most numerical techniques, calculating the Jacobian (most probably by numerical finite difference estimation) is part of solving the equations set up by the harmonic balance method [START_REF] Gasch | Strukturdynamik II[END_REF][START_REF] Genta | Vibration of Structures and Machines[END_REF]. It will be shown that the availability of the Jacobian means that the HBM lends itself nicely to studies of the stability of a solution in the frequency domain, i.e. without having to go back into the time domain (Floquet analysis). To the authors' knowledge, such an algorithm for a non-linear system is novel, as only examples of applications to linear time-variant systems have been found in the literature. Naturally, the algorithm described here is not only suitable for rotor/stator contact problems but also applicable to other common non-linear elements in structural dynamics.

HARMONIC BALANCE FORMULATION

Given the computing resources, the HBM is easily applicable to problems with a large number of degrees-of-freedom (DOFs). Typically, such a problem consists of finite element models of large parts of the structure or substructures where a linear representation is adequate, and some 'problematic' DOFs for special areas, for example where friction, impacts, or other interaction occurs. Usually, the linear DOFs outnumber the nonlinear ones by a large ratio.

The example that will be used later on deals with the dynamics when rotor and stator come into contact. At this stage it does not matter which of the rotor/stator contact scenarios described in the introduction is under investigation. The rotor and stator are modelled as linear structures, and there will be some linear external forces like gravity and out-of-balance. The contact region supplies the nonlinear forcing at a few degrees of freedom on both rotor and stator.

For simplicity, the complete system is split into its linear part, represented by the usual mass, stiffness, damping matrices, with some linear external force vector, f u (for example unbalanced masses), and its nonlinear part, which is represented here as a single force vector f c combining all nonlinear effects (contact between rotor and stator):

[M ] r + [C] ṙ + [K] r = f u (t) + f c (r) (1)
The harmonic balance method offers an alternative to time-domain methods for analysis of cases where a steady-state, periodic solution to the equation of motion is sought. The general idea is to represent each time history, r(t), by its frequency content, R(ω), to obtain a set of equations by balancing the terms with the same frequency components and to start an iterative procedure to find the roots of these equations.

An integer variable, ν, is introduced to accommodate possible sub-harmonics of an external excitation frequency, Ω (for example shaft rotation). The displacements, r(t), and forces, f (t), are represented as truncated Fourier series with N harmonics:

r(t) = N n=1 R n e i nΩ ν t f c (t)= N n=1 F c n e i nΩ ν t f u (t) = N n=1 F u n e i nΩ ν t (2)
Substituting these expressions, (2), into the rotor equation of motion, (1), and balancing the harmonic terms yields, for a harmonic n:

- nΩ ν 2 [M ] + i nΩ ν [C] + [K] R n = F c n + F u n (3)
Bringing all N harmonics into one equation can be symbolised as

K R -F c -F u = 0 ( 4
)
where K is a complex block-diagonal matrix of the following form,

- nΩ ν 2 [M ] + i nΩ ν [C] + [K] = Kn,n (5)
and R and F are the vectors of Fourier coefficients of displacements and forces, respectively. As the Fourier coefficients, F c n , of the non-linear forces, f c , are functions of the displacements (and thus their respective Fourier coefficients),

F c n = F c n (R 0 (ω 0 ), R 1 (ω 1 ), . . . , R N (ω N )) (6)
equation ( 4) is non-linear and must be solved iteratively. This iteration process [6] can be sketched as:

R(ω) (k) FFT -1 ----→ r(t) (k) → f c (t) (k+1) FFT ---→ F c (ω) (k+1) → R(ω) (k+1)
The fast Fourier transform (FFT) and its inverse could be replaced with a (costlier) curve-fitting algorithm that does not rely on an integer number of cycles so that the procedure also works with incommensurate frequency components (quasi-periodic solution). The roots of equation (4) are found using a quasi-newton algorithm, variants of which are widely available in Fortran libraries on the internet or in Matlab's optimisation toolbox. Finite element models of rotor/stator structures can contain quite a large number of degrees-of-freedom. Setting up equation (4) then leads to a much bigger problem with 2N +1 times more unknowns (real and imaginary components for N harmonics and a DC component). Any reduction of the original problem, therefore, leads to a huge saving in computational cost.

REDUCTION

The harmonic balance method offers an elegant means of reducing the problem order, so that only the non-linear DOFs need to be kept [6]. Clearly being able to keep only the non-linear DOFs vastly increases the speed in cases of linear structures which have a few additional non-linear elements, as is typical for many classes of problems. The equation ( 3) is re-ordered for every harmonic, n (the subscripts n are omitted in this section for clarity):

Kmm Kms Ksm Kss R m R s = F c 0 + F u m F u s (7)
where subscripts m and s stand for non-linear (master) and linear (slave) degreesof-freedom, respectively. It follows that

ˆ K = Kmm -Kms Kss -1 Ksm (8) ˆ F u = F u m -Kms Kss -1 F u s (9) and ˆ K R m -F c -ˆ F u = 0 (10)
In contrast to the widely-used Guyan reduction, equation ( 10) is an exact reduction of the original problem as long as the prerequisites for applying the harmonic balance method are met and the number of harmonics included in the decomposition is sufficient. This reduction makes the repetitive task of calculation solutions at various parameters for the continuation of a solution branch much cheaper indeed.

ARC-LENGTH CONTINUATION

Usually, the system behaviour is of interest over a range of values for at least one parameter (for example speed of shaft rotation), so that the solution has to be calculated at different parameter values consecutively. As time-domain methods seem to be used more often than frequency-domain methods, the use of continuation schemes in conjunction with the harmonic balance method appears to be not as common in engineering applications as, say, with shooting methods. However, arclength continuation is just as applicable in the frequency domain as it is in the timedomain. The task of finding a periodic solution for equation ( 1) can be transformed into an equivalent root-finding problem, for example by means of finite difference, shooting, HBM (see equation ( 4))

F (y, Ω) = 0 (11) where Ω is an independent, externally controlled (exogenous), scalar parameter, such as excitation frequency or rotor speed, over a range in which the solution to equation ( 1) is of interest.

Using the notation F y = ∂F (y,Ω) ∂y and F Ω = ∂F (y,Ω) ∂Ω , the differential of equation ( 11) can be written as:

dy dΩ = -(F y ) -1 F Ω (12)
Choosing Ω directly as a continuation parameter (called "sequential" or "natural continuation") fails at turning points because of the singularity of F y (see transition from stable to unstable at Ω ≈ 1.4 in Figure 2). To overcome this limitation, a continuation parameter, α, along the arc-length of a solution branch is chosen, so that y = y(α) and Ω = Ω(α). From equation ( 11) one obtains:

F y dy dα + F Ω dΩ dα = F y F Ω y Ω = 0 (13)
The arc-length, α, may be normalised [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF], so that the tangent vector, y Ω , has unit length:

dy 1 dα 2 + • • • + dy N dα 2 + dΩ dα 2 = 1 (14)
With y(α [k] ), Ω(α [k] ) denoting a solution previously calculated during continuation, an additional equation can be gained for the additional unknown α by multiplying equation ( 14) with (dα) 2 , 0 = g(y, Ω, α)

= y 1 -y 1 (α [k] ) 2 + • • • + y N -y N (α [k] ) 2 + Ω -Ω(α [k] ) 2 -α -α [k] 2 (15)
Equation ( 15) can now be used to parametrise equation ( 4) by formulating an extended system [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]:

F (y, Ω, α) = F (y, Ω) g(y, Ω, α) = 0 (16)
Equation ( 16) allows the solver to move along the arc-length α of a solution branch. The reader is referred to [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF][START_REF] Parker | Practical Numerical Algorithms for Chaotic Systems[END_REF][START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF] for predictor-corrector and step control schemes that facilitate the following of a solution of F (y, Ω, α) along a range of values for α. For simplicity, the results of the previous step y(α [k] ), Ω(α [k] ) may be used as an initial guess for the next, and a simple step control α [k+1] = α [k] + p q can be employed, where p is a constant and q is the number of iterations needed for the previous solution. With quasi-Newton solvers an optimal ratio of p q is suggested by Seydel [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] such that roughly q = 6 iterations of the quasi-Newton solver are necessary in between arc-length steps. As a first guess initiating the continuation procedure, one may use a random guess or in more diffcult cases find approximate values from a solution obtained previously by time integration. As demonstrated in the numerical example below, there are now no problems in passing turning points and following the overhung part of the solution branch shown in Figure 2. Obtaining information about the stability of such a solution branch being followed is discussed in the following section.

STABILITY

When HBM with arc-length continuation is employed, there is nothing in the algorithm per se that can warn the user that a particular solution branch followed has stepped over a turning or bifurcation point and the solution has switched stability, from stable to unstable, or vice versa. For example, there might only be a little change in the conditioning of the Jacobian of the system before and after such a change, nor is there a change in the convergence behaviour of the algorithm. This is a practical problem (not a theoretical one, as the Jacobian F y is indeed singular exactly on a turning or bifurcation point) as change in conditioning of the Jacobian could occur far more rapidly than the step-length is able to resolve.

Stability in the time-domain is usually determined by the well-known Floquet multipliers, which are extensively covered in the previously cited text books. It makes sense to employ this method when one operates in the time-domain anyway. For example by using the shooting method, one can determine the stability in parallel to solving for periodic solutions, albeit at the cost of doubling the problem size [START_REF] Sundararajan | Transactions of the ASME[END_REF]. However, it is less advantageous when one operates in the frequency domain using the harmonic balance method. Many papers employing HBM either ignore the stability issue or revert back to time marching solutions to answer this question. At little computational cost, however, stability can be analysed in the frequency domain with a modification of an algorithm employed for linear time-variant systems. The algorithm is called Hill's method and transforms a linear time-variant system into an eigenvalue problem of a linear time-invariant system [START_REF] Gasch | Strukturdynamik II[END_REF].

In order to use the same approach for non-linear systems, the stability analysis is carried out by investigating the effect of a perturbation around a periodic solution r (t). Although Hill's method for linear systems is documented in text books and the approach to study the stability of nonlinear systems typically involves the linearisation around that solution [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF], the following approach of the combination of the two has not been previously encountered in the literature. Let the perturbation be described as p(t), where p(t) consists of a decay term e λt and a periodic term s(t) [START_REF] Gasch | Strukturdynamik II[END_REF][START_REF] Genta | Vibration of Structures and Machines[END_REF]:

p(t) = e λt s(t) s(t)= N n=-N S n e i nΩ ν t (17) Substituting r(t) = r (t) + p(t) ( 18
)
into equation ( 1), one obtains

[M ] ¨ r + [C] ˙ r + [K] r + λ 2 [M ] s + λ 2 [M ] ˙ s + [C] s + [M ] ¨ s + [C] ˙ s + [K] s e λt = f u (t) + f c r + s e λt . ( 19
)
By substituting the Fourier representations of r (t), s(t) into equation ( 19), the harmonic components can be balanced in an analogue manner to equation (4):

K R + λ 2 M + λ C + K S e λt = F u + F c R + S e λt ( 20
)
where M , C are constructed in a similar manner to K in equation ( 5) and R (ω) , S(ω) are the vectors of Fourier coefficients for r (t), s(t), respectively. In what follows, an attempt is being made to find a cost-effective linearisation for the term F c R + S e λt so that equation (20) can be developed further. Consider a variant of equation ( 4):

F c = K R -F u -E(R) ( 21
)
where E(R) is the error in the balancing terms. Developing this as a Taylor series around a known solution of equation ( 4), R , one obtains

F c (R) = K R -F u -E (R ) R -R + higher order terms (22)
with the abbreviation

E = ∂E ∂R .
Substituting equation ( 22) into (20) and neglecting terms of higher order, equation (20) simplifies to the following eigenvalue problem:

λ 2 M + λ C + E (R ) S = 0. ( 23
)
It is important to note that the term [E (R )] is already available as a by-product of the quasi-newton solution technique, as it is the Jacobian of the objective function defined by equation (4). In the example below, it was sufficient to approximate the Jacobian numerically by calculating the finite differences in the objective function due to small perturbations (a standard implementation as part of a quasi-newton solver). One is, of course, free to use other methods of determining the Jacobian, ill-conditioned cases might require more sophisticated (robust) methods in order to aid the solver to find the periodic solutions (which is a more general issue of solution convergence than HBM or stability analysis).

Solving for the eigenvalues of equation ( 23), one obtains a set of λ i with real and imaginary parts, where a negative real part indicates stability of the solution, as the perturbation p(t) decays with time, and a positive real part indicates instability. So by solving this eigenvalue problem at the end of the overall iteration procedure, and simply checking if any λ i possesses a real part > 0, one can easily determine whether a periodic solution r (t) is unstable. This also helps with finding possible bifurcation points. A change in stability of a solution branch is a sufficient indicator that a turning or bifurcation point has been passed, and the algorithm could be directed to determine the cross-over point within this interval of change more closely. Should this be of interest and the cross-over point found, the rank of F y and F y F Ω at the cross-over determines whether the point in question is a turning or bifurcation point [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. If indeed it is a bifurcation point a further solution branch may be followed.

NUMERICAL EXAMPLE

Although the algorithm described above is easily scalable to systems with a large number of degrees of freedom, a simple modified Jeffcott rotor (Figure 1) is used here for clarity. The equations of motion for a Jeffcott rotor interacting with a linear stator structure are:

0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
m r rr + c r ṙr + k r r r = -f c + Ω 2 m r m e iΩt (24) m s rs + c s ṙs + k s r s = f c (25)
where r r , r s are the rotor and stator displacements in the complex plane and f c is the contact force between rotor and stator,

f c = k c δ ( 26
)
where δ is the depth of the contact described below in (27) and k c the local (in this case linear) contact stiffness (for more realistic simulations one would have to choose a nonlinear contact force, see [START_REF] Fumagalli | Sixth International Conference on Vibrations in Rotating Machinery[END_REF]). For the purpose of numerical simulation, a small contact penetration δ of the rotor and stator rings is allowed. The contact stiffness k c in this penetration region is being set to a value orders of magnitude higher than the rotor shaft or stator suspension stiffness, so that the penetration depth is orders of magnitude lower than rotor and stator deflections. The contact depth is defined as

δ(t) = r r + r e iΩt -r s -s -h e iψ if |r r + r e iΩt -r s -s | > h, 0 o t h e r w i s e ( 27
)
where r r , r s are rotor and stator displacements, h is the gap size, r a possible offset of the rotor disc and s a stator offset. These entities are depicted in Figure 1.

For the special case of full annular rub with r , s = 0 and isotropic rotor supports, the equations of motion become quasi-static for pure forward or backward whirl. At a given speed, the steady-state conditions of rotor whirl are such that the radial deflection of the rotor is constant. The only frequency component in the imbalance response spectrum is thus the engine-order speed, Ω, thus N, ν = 1. This simple case is used here to illustrate the stability and continuation study, with the following values: m r = 1, k r = 1, c r = 0.02, m s = 0.1m r , k s = 2k r , c s = 0.002, m = 1, friction μ = 0.1, contact stiffness k c = 100k r (1 + iμ), gap h = 3. In Figure 2 the magnitudes of the rotor and stator responses r r , r s are plotted versus the rotor speed of rotation Ω. One can see that at speeds Ω < 0.9 the rotor imbalance response is too low to overcome the clearance (h = 3, dashed line) and rotor and stator are not in contact (stator response zero). As the continuation procedure is started in this speed range with no rotor/stator contact, a simple random guess for the initial conditions is sufficient for solver convergence. At speeds 0.9 < Ω < 1.4 rotor and stator are in contact (non-zero stator response), albeit the overhung part of the curve represents an unstable solution. At speeds Ω > 1.4, well past the natural frequency of the rotor, which has been normalised to ω r = 1, the super-critically running rotor loses contact with the stator. Figure 3 shows a second solution branch at Ω > 6, which is not seen in the 1 DOF Duffing type oscillators that display only the overhung behaviour in Figure 2. It must be noted that by following the branch previously discussed, the one that lost contact with the stator and is coming into this picture from the left, there is no indication of the existence of the second branch. The branches of the rotor deflections intersect, but in the whole space of rotor and stator deflections these curves do not come near each other (see stator solution branches), so there is no warning in terms of changes in stability or conditioning of the Jacobian.

This second solution branch was found by brute force, using hundreds of random initial guesses at different speeds, in a quest to find out whether the system could vibrate in an 'inverted' modeshape, where the stator is moving like a hola-hoop around the rotor. The success rate of finding the second branch from the random initial guesses as a whole were poor: many did not converge to any solution, most settled on the solution where rotor and stator were out of contact. The ones that did settle on the second solution branch seemed to settle just as easily in terms of convergence on the unstable part of the branch as on the stable one, albeit a physical system would object to that. Once a solution on that branch is found, it is easy to follow with a continuation scheme.

Making the system more general with r , s = 0 and non-isotropic rotor and stator mounts gives rise to more complicated motion. Under these circumstance the system has speed regimes where the rotor is rubbing along the stator only over parts of the orbit, causing sub-and super-harmonics to emerge. Even a full annular rub in these circumstances would consist of higher harmonics (but no sub-harmonics, which require intermittent contact). Figure 4 shows the frequency content of the rotor motion at various speeds. The harmonics are expressed in engine order ratios and their magnitude is given in dB. At lower speeds (<15Hz) only the 1EO (= first engine order = Ω) component is present, as rotor and stator are out of contact and without the nonlinear forces no other frequency components but 1EO (due to the imbalance) is expected. As soon as rotor and stator come into contact at higher speeds do the other frequency components start to come into play. ν = 12 and n = 4, 6, 12, 24, 36, 48, 60 are chosen to capture the frequency components 1/3, 1/2, 1, 2, 3, 4, 5 EO.

These frequency components are easily confirmed by comparing the results with

m r = 1, k r = (2π17) 2 , c r = 4, m s = 2, k s = (2π12) 2 , c s = 3, k c = 100k r , h = 2, m = 0.2, r = 0, s = 1.
). There is no visible difference between the orbits obtained with time-marching and HBM methods, so for the purpose of this example one may regard the HBM solution with 32 harmonics as exact. However, 32 frequency components is a large number (which are chosen here only to validate that there can be complete agreement between timemarching and HBM results), and it is interesting to see what happens when fewer frequency components are used. Figure 5(b) shows the resulting orbits when the HBM is used with only 4 frequency components in the setup (n = 1, . . . , 4, ν = 2). Not surprisingly, the orbits have a slightly different shape compared to Figure 5(a), but one can see that the fundamental physical effect, namely the sub-harmonic loop, is still retained. When reducing the number of harmonics further to two or only one, the solver only found a physically totally different solution, namely rotor and stator out of contact, where the stator is stationary and the rotor in a simple out-of-balance orbit. This emphasises the problem to set up the equations to contain those frequencies that the system responds at. It is therefore a method that in general can only be used if some a priory knowledge about a system is available, be it from experience with previous designs, time-marching simulations, or experiments. Measurements from a test rig [START_REF] Ewins | Seventh International Conference on Vibrations in Rotating Machinery[END_REF] also show that the physical system in such a rotor/stator contact scenario predominantly responds at a few distinct frequencies, thereby indicating that the HBM with a few selected frequencies is an adequate representation of the physical behaviour.

CONCLUSIONS

The presented frequency-domain algorithm calculates the periodic solution, traces the solution along a varying parameter and determines the stability of a solution branch. For the investigation of stability previous studies have reverted back to the time-domain. A major motivation for working in the frequency domain is the computational speed advantages it has over time-domain methods.

The solution process of the HBM method itself is less expensive than time-domain methods, and the reduction to the non-linear degrees-of-freedom offers vast savings in large finite element models with only a few non-linear components. It was found that in the given numerical example the harmonic balance method was over 100 times faster than time-domain shooting and boundary-value-problem solving.

HBM also functioned properly in many instances where the time integration routines had difficulties. The source of these difficulties lies in the contact problem. As both rotor and stator have non-negligible mass (and thus dynamics in their own right), the penalty stiffness k c determines the violence of the rotor/stator impacts and thus tests the robustness and influences the speed of convergence of any chosen method. A more sophisticated contact model can alleviate this problem to a large extent.

The above mentioned advantages of the algorithm used here make it a worthwhile ingredient in any tool set analysing non-linear systems. However, there are, of course, situations where the harmonic balance method does not find a solution. This is usually an indication that not enough or not the right harmonic components are included in the setup. Furthermore, the fact that a solution has been obtained is not a sufficient condition that no important frequency components have been omitted in the setup. One future enhancement might be to choose a hybrid approach of an initial setup in time-domain, giving an overview of the response spectrum, and continuation in frequency-domain to maximise on the robustness and speed advantages of the HBM.

Figure 2 .

 2 Figure 2. Magnitude of rotor and stator response at constant speed. solution: + stable, o unstable

Figure 3 .

 3 Figure 3. Second solution branch at Ω > 6. solution: + stable, o unstable

Figure 4 .

 4 Figure 4. harmonics 1/3, 1/2, 1, 2, 3, 4, 5 EO versus speed of rotation [Hz]

Figure 5 .

 5 Figure 5. HBM rotor and stator results with various numbers of harmonics in the setup

ACKNOWLEDGEMENTS

The authors are grateful to Rolls-Royce plc. for providing financial and technical support for this project and for giving permission to publish this work.