Contacts in multibody systems

Research in contact mechanics is reviewed, mainly based on result determined at the Institute B of Mechanics in the Thchnical University of Munich. Interest in such problems arise in the 1960s in relation to statics and elastomechanics. Most of the mathematical tools were developed and used in those fields. As a result of increasing pressure from the practical side concerning vibrations and noise related to contract processes, methods of multibody dynamics with unilateral constraints have been elaborated.


INTRODUCTION
Contact phenomena, including bilateral and unilateral ones, often arise in dynamical systems. Walking, grasping and climbing are typically unilateral processes; the operation of machines and mechanisms includes a large variety of unilateral behaviour. From this the need arises to extend multibody theory by analysing contact phenomena.
All contact processes have some characteristic features in common. If a contact is closed, the motion changes from slip to stick, which leads to some additional constraints generating constraint forces. We then call the contact active, otherwise it is passive. Obviously transitions in such contacts depend on the dynamics of the system under consideration. The beginning of such a contact is indicated by kinematical quantities like relative distances or relative velocities and the end indicated by kinetic quantities like the change in normal forces or friction forces. This serves as a basis for the following mathematical formulations.
A large variety of possibilities exists in modelling local contact physics, from Newton's Poisson's and Coulomb's laws to the discretization of the local behaviour by the FE-or BE-methods. However, simulations of large dynamical systems require compact contact laws. Therefore, we shall concentrate on the first type of laws which, despite their simple structure, are still able to describe a large range of applications realistically.
The literature covers aspects like contact laws, FEM and BEM analysis contact statics, contact dynamics and a large body of various applications, (see for example, [1]). As regards multibody systems with multiple unilateral contacts, most of the mathematical fundamentals, though initially regarding static problems only, were laid down by European scientists. The first investigations were carried out by Moreau and his school (2] and future time been developed considerably [3]. Moreau introduced convex analysis into multibody dynamics and reformulated the classical equations of motion in terms of measure differential inclusions in order to cover both, impact free motion and shocks, as they occur in frictional contact problems.
Powerful methods for static and dynamical problems of contact mechanics were developed in (4,5]. Panagiotopoulos in particular established a general theory of unilateral problems in mechanics (4]. LOtstedt developed an advanced index-two-type integration algorithm for planar contact problems in rigid-body systems (5]. The Swedish school in this field has continued with remarkable results by Klarbring, who focuses his work on problems of FEM-and BEM-modelling [ 6]. Russian scientists have developed [7,8] methods for non-smooth transformations of variables, which provide regularization in systems with impacts.
At the Institute B of Mechanics in the Technical University of Munich research has been carried out in this field for more than ten years, which is mostly summarized in [9]. Some more recent results on nonlinear complementarity problems may be found in [10]. As the Institute B of Mechanics is one of engineering mechanics, most of the research work deals with transferring the rigorous mathematical fundamentals to an engineering and application-friendly level. 2

THE EVOLUTION OF THE THEORY
In the following we present some cornerstones in the evolution of dynamic contact theory.
1. The woodpecker toy. The first example with structure-variant properties, which was analysed at the institute, was a woodpecker toy which operates by self-excited vibrations [11 ]. Figure 1 shows a mechanical model of the toy and its operating principle. Here, mm and lm are the mass and central moment of inertia of the sleeve, and ms and ls are the corresponding quantities of the woodpecker. The angles between the vertical and the axes of the sleeve and the woodpecker are denoted by 'PM and 'Ps· They take values 'Pkt and 'Pk2, respectively, when the sleeve or the beak touches the pole. At these values the constraints are switched on and off.
Self-exited oscillations of the woodpecker are caused by the presence of some energy source ST, here gravitation. In the absence of contact events, they are described by the functions included in VI. The contact events themselves act as switching conditions initiated by the sleeve as the main switch SS, and the beak as an auxiliary switch SB, not necessarily needed for the toy to work. These switches control the energy transfer from gravitational energy as the energy source ES into kinetic energy in the form of linear and angular motion, as well as into strain energy stored in the spring. The terms ri denote different types of feedback acting in the system, such as the switching on and off of unilateral constraints (rt> r 2 ), energy transfer by impacts (r 3 ), and the dynamical couplings between the coordinates r 4 • The theory employed in 1984 was at that time incomplete because no reasonable concept for impacts with friction was available. Friction was therefore included empirically. The results a great very well with measurements.
In subsequently an impact-with-friction theory appeared [12,13]. Its application to the woodpecker example confirms the simpler theory while at the same time improving the friction model. Figure 2 shows some results in the form of phase portraits [12]. The data used for the woodpecker may be attained    ... from [9]. The main phases of the motion are ( Fig. 2): 1-2 sliding, 2loss of contact, 2-3free fall with high-frequency oscillation; (73Hz), 3inelastic impact of the upper edge of the sleeve, 4beak impact, 5-6free fall downward, 6-7transition to sleeve sticking and self-locking, 7-1woodpecker angular motion with a low-frequency oscillation of around 9 Hz, which has also been confirmed by measurements.
2. Gear rattling. Gear rattling is a noise problem which occurs in gear trains not under load [14].
Examples are change-over gears or gear-driven balancers. Figure 3 illustrates a five-stage gear configuration where the drive shaft AS moves the countershaft CS, and according to the switched gear, the countershaft drives the main shaft ES. All gear wheels mesh. Switching is performed by synchronizers connecting the appropriate wheel with the main shaft. The wheel is loaded and all other wheels are not under load and may rattle due to backlash in the mesh of the gears. The drive shaft is excited by the torsional vibrations of the drive system, which are more or less harmonic. These excitations are transferred to the countershaft where they produce rattling in the loose gear wheels. Figure 3 shows the case with the fourth stage switched, which for this type of gearbox results in a direct connection of the drive and main shafts. The countershaft and all gear wheels are not under load. The harmonic excitation leads to an impact driven vibration of the countershaft (below on the left) and, for example, in the fifth gear, to chaotic vibration of the fifth gear wheel (below on the right) [15]. An excellent relative measure for the noise generated by gear rattling can be achieved by summing up all impact impulses AN and introducing the noise intensity CL I AN I· Many applications to real gears  confirm that this measure describes all parameter dependence correctly. Figure 4 shows two examples for the change-over-gears of a pick-up truck (1) and a higher category medium-class car (2). Only one measurement point has been considered for choosing the proportionality factor c between the measured noise intensity (solid curves) and the calculated noise intensity (dotted curves).

Turbine blade damper.
The two examples of the woodpecker toy and gear rattling only demonstrate impact theory, supplemented by friction. The operation of a turbine blade damper is entirely based on stick-slip processes, which was at that time a step in a new research direction [16,17].
Thrbine blade dampers are parabolic sheet steel devices, often used in gas turbines to damp the blade vibrations by relative motion through dry friction. Figure 5 shows a typical configuration. It consists of two platforms 2 with masses m 1 and m 3 , each carrying a turbine blade 1 with mass m 2 and m 4 , respectively.
The platforms are separated by a frictional damper 3 of mass m 5 and moment of inertia J 5 • The displacements of the platforms and the turbine blades are described by the coordinates qt> q 2 , q 3 and q 4 • The unconstrained damper has two translational degrees of freedom q 6 and q 7 and a rotational one q 5 • The constant radial force Fz represents the centrifugal force and presses the dampers against the oblique planes of the platforms. The blades are excited approximately harmonically by gas-dynamic forces F 2 (t) and F 4 (t). The elastic and dissipative properties are modelled by four pairs of springs and dashpots ci and di.
For such a model the following motions are possible: sliding at both contact points K 1 and K 2 (5 DOF); sliding at K 1 and stiction at K 2 ( 4 DOF); stiction at K 1 and sliding at K 2 ( 4 DOF); stiction at K 1 and K 2 (3 DOF). Stiction means for the 4 DOF case rolling without sliding, for the 3 DOF case a reduction of the coordinates (qt> q 3 , q 5 ) to 1 DOR For each configuration a set of equations of motion can be established, which describe the motion as long as no transition occurs. In the case of a blade damper these transitions are typical of stick-slip behaviour.
Transitions from sliding to stiction occur when the relative tangential velocity at one of the points of contact becomes zero, and at the same time the static friction force must be larger than the tangential constraint force.
Transition from stiction to sliding occurs when the tangential constraint force becomes equal to the static friction force, and at the same time the tangential acceleration is non-zero.
After each transition event the appropriate set of equations of motion must be selected to describe the further behaviour of the system up to the next transition [16]. Figure 6 shows the dependence of the excitation force F 24 (in Newtons) on the frequency [ 24 (in Hertz) for three different values of the centrifugal force Fz and for a platform angle 'Y = 50°. The amplitude of the exciting force F 24 must be increased for larger centrifugal forces Fz. The very small F 24 values for an excitation frequency of 61Hz can be explained in the following way. The value of/ 24 =61Hz corresponds to the natural frequency of the blocked-damper system. This means that an excitation at this frequency leads to large oscillation amplitudes of the blades and as a consequence to a considerable reduction of stiction. Therefore, the damper devices move from stick to slip when excited by a force of small amplitude.

4.
Landing impact of an aeroplane. The landing impact of an elastic aeroplane is governed by the configuration of the aeroplane and its elasticity and by the main and nose under carriages. The practical example of a landing aeroplane was the first one where the different sets of equations of motion could not be written down explicitly, at least without considerable effort, but where some theoretical algorithm had to be established which enable the appropriate equations to be shown automatically [18,19]. The internal structure of the undercarriage systems leads to impacts and to stick-slip phenomena, which influence the landing dynamics significantly. A plane elastic model was established [9] for the following bodies (see Fig. 7): (1) and elastic fuselage (Sis its centre of mass), (2) an elastic main undercarriage, (3) an elastic nose undercarriage, (4) rigid main undercarriage wheels, (5) rigid nose undercarriage wheels. For the shock absorbers it will be assumed that their bending deformation is homogeneous over the variable length, which means that the two components of the absorber do not bend in a different way. This turns out to be a sufficient approximation.
During operation we have to consider various continuous and discontinuous non-linearities. They tyres possess a stiffness characteristic which increases progressively with the tyre deformation. Moreover, the spin-up process of the wheels is governed by a highly non-linear friction-slippage relationship. The shock absorber itself reveals some non-linear features. The two gas chambers follow a polytropic compression and expansion law, the fluid flow losses induce a quadratic damping behaviour, and the friction forces from the ground reaction forces are non-linear also.
Discontinuous non-linearities arise but to the presence of upper and lower stops within the shock absorber and through stick-slip processes between the two shock absorber cylinders. The various discontinuities in the overall system would have required more than 40 different sets of equations of motion, which was not feasible. At present an algorithm for the systematic and automatic generation of these equations has been developed, based mainly on the properties of the constraint matrices W (see Section 3 below and [9]). The results agree well with practical experience.

Assembly processes.
Mating different parts together involves various and sometimes quickly changing combinations of impacts with and without friction, of stick-slip processes and even of jamming. As these combinations and the transitions between them cannot be predicted in advance, because they depend on the current state of the system, we need a theory which can deal with such processes. It was this problem that led to the development of the methods of complementarity, sub-differentials and convex analysis [9,12,20].
The practical application of these methods to some extent, abstract theories usually operates on different levels: the transition between contact configurations is either indicated by magnitudes of relative kinematics (at the beginning) or by those of constraint forces (at the end). The relative kinematical magnitudes undergo a metamorphosis from indicators to constraint equations as components of the constraint matrices, valid for the new contact configuration. The transition itself is governed by the principle of complementarity of the kinematic characteristics constraint forces: if one group vanishes, the other does not. Considerations of this kind in combination with the accompanying mathematics [4,9] allow of efficient modelling of the assembly processes. As an example we take the operation of inserting a peg into a hole, performed by a manipulator [21,22].
A robot manipulator is usually considered as a tree-like multibody system with rigid and elastic components. Links are very often modelled as rigid bodies and only for special tasks as elastic bodies. Joints may ideally follow a given torque law, or they may be elastic with their own degrees of freedom (see Fig. 8). Therefore, we may subdivide our total number n of degrees of freedom into motor or internal degrees of freedom (nM), external ones (nA), and elastic ones (nc)· The equations of motion are then derived by applying Jourdain's principle which says that constraint forces at the joints are perpendicular to the direction of motion and do not affect the power.
Because most mating tasks in assembly cells are limited to a small area of workpiece interaction, the robot motion will be slow and centrifugal and Coriolis forces can henceforth be neglected compared to gravitational and inertia forces. Hence, the robot dynamics can be linearized around a given system state. During assembly the workpiece held by the robot comes into various contact situations with the complementary mating part. This is illustrated in Fig. 8 for a planar peg-in-hole insertion task. Here the following notation is used: 1 is the arm, 2 is the motor, 3 is the gear model, 4 is the compliance element, and 5 is the gripper. In this process one or more contacts may arise which constrains the gripper motion and thus the robot's mobility by closing the kinematic loop of the manipulator.
Experimental investigations have focused on a verifying of the theory for various configurations of assembly processes and on a considering perturbations during the peg-in-hole insertion. We will first Fig.8 consider the perturbed motion during a planar peg-in-hole insertion task where forces only develop in the presence of uncertainties. To verify our theoretical approach, experiments with a five-degree-offreedom laboratory robot with a force-torque sensor were carried out. The fixture housing the complementary part is equipped with six distance sensors that measure the gripper's position and orientation and in addition disposes of two translational and three rotational adjustments which enable definite positioning errors to be produced with respect to the parts that are to be mated.
In one experiment we investigated the insertion process when there was a lateral error between the peg and the hole. The respective numerical results (on the right) and experimental results (on the left) are shown in Fig. 9, where t is time and Fx is the insertion force during assembly due to lateral offset.
The peaks in the force history at t = 0.9 s result from the impact when the peg hits the chamfer. When two-point-contact occurs at t = 1.2 s, the manipulator's motion is slowed down very rapidly until a jamming condition is reached. Only when the controller torques become large enough, does the motion continue with sliding at two and finally even at one contact point. The cycle represented by the oscillations in the force history repeats until the end of the trajectory and demonstrates the phenomenon of changing constraints during an assembly operation. 3. MATHEMATICAL FORMULATION Consider a system of rigid bodies with nA frictional contacts. The maximum number of degrees of freedom f is obtained when none of the contacts is closed. In this state the system may be described by a set of generalized coordinates q E Rf, and the equations of motion take the form Here, M is a symmetric and positive-definite mass matrix, h contains the gyroscopic: accelerations and the single-valued applied forces, and the sum relates to the contact forces. For each pair of contact points i we take an orthogonal system (n, u, v) such that n becomes the common normal to the surfaces of the contacting bodies and (u, v) specify two orthogonal directions in the tangent contact plane. The triple (n, u, v) then defines three connections for the generalized forces wN and Wr = (wu, wv) in the configuration space, and "-N together with "-f = ("-u, "-v) are underline to be the scalar values of the normal and the two tangential projections of the contact forces, respectively.
The kinematic state of a contact is determined by the distance 8N(q, t) between the contacting, bodies by the relative velocities (gN, gu, gv) of the contact points in the three directions (n, u, v) and by their time derivatives (jjN, gu, gy). A straightforward evaluation of the contact kinematics yields (19,21 ]. (3.2) with the same WN, Wras used in (3.1) In (3.2) 8N is the distance between the contacting bodies, gN describes the relative velocity of the contact points in the normal direction and gT = (Eu, gv)T is the tangential relative velocity, which also determines the sliding direction of the contact points in the common tangent plane relative to each other. Differentiating 8 with wN = wf4li-N and wr = w[q + ~-This equation is needed to formulate the frictional contact laws on the acceleration level, in the same way as when reducing the index of a differential algebraic mechanical systems from 3 to 1. Generally, every force occurring in a rigid multibody system may be expressed in terms of some relative kinematic magnitudes such as the relative displacements and velocities. Classical applied forces are represented by continuous functions whereas classical bilateral constraints already require a representation via set-valued mappings. Contact laws, as used below, belong to the latter class and constitute mixed-type force characteristics, acting in some areas as classical applied forces and in other areas as constraints.
Suppose gN; is the distance between the contact points and liN; is the corresponding normal force.
We may then express the non-deformability of rigid bodies by the Signorini-Fichera: condition [3][4][5][6]12] (3.4) which can be combined, for example, with Coulomb friction in the tangential directions [2][3][4]6], i.e. By using certain continuity assumptions on the trajectories q(t) can also express the contact laws (3.4), (3.5) on the acceleration level in order finally to solve Eqs (3.1) for the unknowns ij. This yields, for the contact law in the normal direction (3.4) The complementarity condition in the second line of (3.7) can also be expressed in the form of one of the following variational inequalities [3,5,6,10]  Similarly we obtain from (3.5) the friction law on the acceleration level in the form [10] Ar. = 0 for i e I A \ IN ' where e; is the unit vector in the sliding direction, i.e. e; = gd I8Ti I· We recall: that, due to the assumed dependence on the normal forces, only quasivariational inequalities are allowable for the third portion of the contact law (3.9), which take the form [10] I.  .7) and (3.9) provide a complete description of impact-free non-smooth rigid body motion under the influence of Coulomb friction within the framework of non-smooth analysis. The equations considered above thus allow of further discussion and evaluation by applying available theoretical results and algorithms from this field [2][3][4]6]. Impacts have been excluded, but they can be treated in an analogous manner by rewriting (3.1) in measure space and solving it at the impact times [2]. This yields, together with appropriate impact laws, a set of relations similar to those considered above, in which the momenta play the role of forces, and the accelerations are expressed by velocity jumps [9,12,13] 4. CONCLUSIONS In the theory of multibody systems with unilateral contacts the methodology has reached a state, which enables to be applied in nearly all fields of mechanical engineering. Since unilateral contacts occur very often in machines and mechanisms, the use of the appropriate theories is pending. This paper gives some examples which demonstrate the typical features of unilateral constraints in combination with machines, e.g. they may ensure certain functional behaviour like the transportation rate in vibratory feeders, or they might generate undesirable vibrations, noise and wear like rattling in gears or wear in chains.
The main difficulties limiting applications today are of a numerical nature. To describe a machine like the vibration conveyor requires tedious numerical evaluations for solving the complementarity problem at each time step. Even with modem high-speed computers severe computing time problem arises. At present there are three groups of algorithms available suitable for calculations.
First, there are methods which try to find a solution by a kind of intelligent trial and error procedure. They are unsatisfactory. The pivot-algorithms are related to the well-known simplex-algorithm. The Lemke method is an example. It can be used with some success, but requires large computing times. Within the framework of iterative methods, a modified Newton approach seems to be very promising, even for the non-linear complementarity problem. It is being developed at the present time.
A second problem, also closely related to numerical calculations consists in the solution of non-linear complementarities, which arise when considering spatial unilateral contacts. Up to now we have used the following methods: linearization of the friction cones, augmented Lagrangian methods and NCPfunctions, which visor in Mangasarian's theorem. All these methods work satisfactorily, but all need improvement, or at least must be developed. We are also working on this. We note, in conclusion that the theories and methods available today help engineers to understand much better problems with unilateral contacts, which formerly were treated only in a very simple way. We are now able to analyse mechanisms with unilateral contacts in a way which, 5 years ago, could not even have been thought of.